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ABSTRACT 


This is an exhaustive study of the unconstrained optimal control of relative motion in 
the orbit of two spacecraft using low thrust. Orbital maneuvers using low thrust is a 
wide field of research with many possible applications. Whereas conventional impul¬ 
sive maneuvers uses thermal rocket engines which produces high thrust acceleration 
for a relatively small period of time, the low-thrust maneuvers can use unconventional 
methods such as electric rocket engines and solar-sail systems which are of relatively 
longer duration and use continuous low-thrust for the operation. In this thesis, mini¬ 
mum fuel optimal relative motion of two spacecraft in low-orbits using Euler-Lagrange 
formulation is studied. Both linearized Hill-Clohessey-Wiltshire model for a spherical 
planet and the model including the approximate zonal harmonics of a non-spherical 
planet are incorporated in this study. The main objective of this study is to examine 
the applicability of the approximate models in the optimal relative transfers. Effects 
of the variations in the initial conditions, that are, initial distance and initial velocity 
of the maneuvering spacecraft relative to the target spacecraft have been studied in 
different cases of the optimal control problem of the relative motion in the orbit. By 
analyzing the results of the numerical solution of the optimal control problem for the 
relative transfer, it has been found that for the same specified time of the optimal 
transfer, as the initial relative velocity of the maneuvering spacecraft increases, the 
error in the optimal control trajectory as obtained using the approximate equations 
by comparing with the non-linear equations, increases. When the relative distance 
between the maneuvering and the target spacecraft increases, it is observed that 
the error in the optimal control trajectory increases . However, keeping the relative 
distance and the initial relative velocities of the maneuvering spacecraft constant, 
and increasing the time of the operation, the error in the optimal trajectory is found 
to be decreasing upt-o a certain time limit and then increases again. By studying 
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various cases of the optimal coplanar stationkeeping problem, it has been found 
that there is a certain limit on the specified time of the relative transfer beyond 
which the maneuvering spacecraft crosses the specified nominal distance relative 
to the target. To check the accuracy of the numerical solution, the closed form 
solution of the optimal relative transfer problem for the spherical Earth model has 
also been obtained and different cases have been studied and the numerical results 
were found to be in good agreement. For studying the effects of variations in the 
initial conditions when oblateness of the planet is also taken into account, different 
cases have been studied in which the comparison has been made with the spherical 
planet model. It has been found that keeping the initial relative distance of the 
maneuvering spacecraft and the specified time of the relative transfer constant, when 
the initial relative velocity of the maneuvering spacecraft increases, the error in the 
control input also increases. The comparison has also been made by varying the 
inclination of the reference orbit, and the results of the numerical study have been 
presented in this report. 
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Chapter 1 


Introduction 


1.1 Relative Motion in the Orbit 

Relative motion in the orbit means the motion of a spacecraft orbiting a planet 
relative to the other spacecraft orbiting the same planet. There can be one primary 
spacecraft known as the target and the other spacecraft with the task of performing 
the required maneuver relative to the target. Based on the mission requirement, 
the various relative orbital transfers can be rendezvous and docking operations, and 
maintaining station relative to the target. 

Rendezvous operation means that the target spacecraft has to perform an orbital 
transfer from its initial position in the orbit such that its relative position and velocity 
relative to the target after the specified operation time should be zero. Docking of a 
maneuvering spacecraft, such as, Soyuz rocket, with the target, such as, International 
Space Station (ISS) is an example of such operation. In the stationkeeping missions, 
it is required for the maneuvering spacecraft to maintain a specified distance from 
the target spacecraft within a specified time-frame. Based up on the type of the 
mission, there can be more than one maneuvering spacecraft also, with each spacecraft 
maintaining the specified relative distance from the target. Such missions are known 
as spacecraft formation flying. As per the mission requirement, it is also possible 
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that the target is virtual and there can be constraints on the relative distances of 
the maneuvering spacecraft, so that, there will be no collisions among them. In a 
spacecraft formation, multiple spacecraft can work together in a formation for various 
purposes such as research, mapping, and communication. A spacecraft formation 
containing multiple small and relatively cheaper set of satellites is an attractive 
alternative for one large costly satellite. 


1.2 Impulsive and Non-Impulsive Transfers 

Conventionally, thermal rockets are used for the various relative transfers in the 
orbit. They uses thrust impulses, that is, large force in a short period of time, to 
instantaneously change the velocity of the spacecraft, pushing it into a new orbit. 
Based on the requirement of the operation, the operation can involve two or three 
impulses. 

Unlike using a thrust-impulse to instantaneously change the velocity of the spacecraft, 
in non-impulsive transfer, there is a continuous application of thrust, so that, the 
spacecraft changes its direction gradually. Non-impulsive transfers relies on the 
low-thrust propulsion for the operation. Presently, there are various methods for 
the low-thrust propulsion, which are also discussed in [2] in detail. Some of the 
mentionable low-thrust propulsion methods are, ionic propulsion, Hall-effect thruster 
and solar-sail systems. The electrostatic ion thruster uses high-voltage electrodes to 
accelerate ions with electrostatic forces, and achieve a specific impulse within the 
range of 4000-8000s [34], Solar-sail systems uses ultra-thin sails in the space which 
gets their motion by the pressure of sunlight on their large surface, without using any 
stored propellent on-board. Solar sails, which are capable of generating a specific 
thrust of the order of 2 x 10 _4 g [2] can offer low-cost operations with a long life of 
the system. 
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1.3 Optimal Control with Low Thrust 

Low-thrust propulsion in the space is emerging as the possible future technology for 
the various space missions, ranging from the intraplanetary orbital operations such as 
rendezvous of the spacecraft with a space station to the interplanetary explorations 
involving reaching to farther planets. Whereas the conventional missions uses 
mostly the thermal rocket engines, which have a considerably high specific thrust 
value and higher mission cost, the low-thrust missions are aimed at using the 
unconventional propulsion techniques such as electric rocket engines and sail systems, 
with a considerably low specific thrust output, thereby, reducing the mission cost to a 
minimum. This thesis is concentrated on the low-thrust orbital missions, with planet 
Earth as the primary body, but the same analyses can be extended to any planet 
with the change in values of certain parameters. Low-thrust orbital missions consists 
of operations like rendezvous and docking operations, station-keeping missions and 
the satellite formation flying missions in the orbit. A low-thrust rendezvous mission 
means that the maneuvering spacecraft i.e., the spacecraft which is using low-thrust 
motion, has to perform a rendezvous to the fixed target spacecraft, such as, the space 
station such that, the final relative distance and the relative velocities should be 
zero. 

The general control system problem can be classified into terminal control systems 
and tracking control systems. In this case, the terminal control refers to the guiding 
of the maneuvering spacecraft from its initial position to the final position relative 
to the target within a specified time of the operation. Once the reference nominal 
trajectory for the states and the control input is obtained, it is required to design 
a tracking controller which can ensure that the maneuvering spacecraft will follow 
the nominal trajectory. The tracking control system is responsible for maintaining 
the required trajectory of the spacecraft as predicted by the terminal control. This 
thesis is focused on optimizing the terminal control. 

Consider the state of the maneuvering spacecraft is defined using the state-space 
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equation as, 


x = Ax(t) + Bu(t) 


( 1 . 1 ) 


in which, u(t ) is the required control input (specific thrust). The objective of the 
terminal control system is to guide the maneuvering spacecraft from initial state 
T(0) to final required state x(tf ) in the specified mission time, tf. The performance 
of the terminal control system for the maneuvering spacecraft is defined using the 
performance index, 


J 


1 

2 


rtf 


u T Ru dx 


’to 


( 1 . 2 ) 


where, R is the positive definite weight matrix. 

For the fuel-optimal control, the above stated performance index needs to be mini¬ 
mized with respect to the control input, u. Since, the state and the control inputs 
are related to each other using the state equation, hence, the optimization must be 
subjected to the equality constraints of the state-space equation. Therefore, the 
augmented function known as Hamiltonian is stated as, 


7-L(x(t),u(t),t ) = u T Ru + A T (Ax(t) + Bu[t )) (1.3) 


where, A is the co-state vector having the same size as the state vector. The necessary 
condition for stationarity states, 


OR 

d u 


= 0 


(1.4) 


Equation (1.4) gives the value of u in terms of the co-state variable A, which is 
then used in eqn. (1.3) to obtain the optimized value of the Hamiltonian, R*. The 
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derivative of the co-state vector is determined as, 



And, the derivative of the state-vector is given by, 



(1.5) 


( 1 . 6 ) 


Now, Eq. (f.5) and Eq. (1.6), which forms the two-point boundary value problem 
can be solved for the given boundary conditions, i.e., x(t = 0) and x(t = tf ) to get 
the solution for the states and the co-states, which can then be used to find the 
optimal control using Eq. (1.4). 

The above procedure produces a reference nominal trajectory of the states, x d and 
the nominal control input trajectory, u d , required for the stated optimality. 


1.3.1 Numerical Solution of the Two-Point Boundary Value 
Problem 

The two-point boundary value problem can be solved by various numerical methods. 
Tewari [16] enlists such methods for solving two-point boundary problems such as 
the collocation method and the shooting method. As discussed in section (1.3.1), the 
two-point boundary value problem for obtaining the optimal trajectory conprises of 
two sets of first-order ordinary differential equations, where the first set is given by 
Eq. (1.5) as the state-equations, and the other set is given by Eq. (1.6) as co-state 
equations. 
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Shooting Method 

Consider the two-point boundary value problem given as, 

x" = f(t,x,x'), to<t<tf (1.7) 

with boundary conditions, 

x(to) = x 0 , x(t f ) = Xf 

The shooting method treats the two-point boundary value problem given above as 
an initial value problem, in which t is the time, having initial and final valules as 
to and tf, respectively. Hence, the above problem can be posed as the initial-value 
problem as, 

x " — f(t, x i x ')i to<t<tf (1.8) 

with initial conditions, 

x (t 0 ) = cco, x> (t 0 ) = y 

in which, y must be chosen in a way that it satisfies the remaining boundary condition, 
x (tf ) — Xf. The shooting method requires the selection of the proper slope so that, 
the resulting solution shoots to x(t) — Xf at t — tf. The ordinary differential equation 
associated with the initial value problem can be solved using standard numerical 
methods such as, Runge-Kutta or multi-step methods (e.g., Adam-Bashforth method 
and Adams-Moulton method). 
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Collocation Method 

Another method used to solve the two-point boundary value problems is the colloca¬ 
tion method [16]. This method approximates the solution by a linear combination of 
a number of piecewise continuous polynomials that are usually defined on a mesh of 
collocation points. Thereafter, the approximate solution is then substituted in to the 
system of the ordinary differential equations such that the system is exactly satisfied 
at each collocation point. The most common choice of approximation is a linear 
combination of spline functions. The use of collocation by splines to numerically 
solve the two-point boundary value problem has been presented by Ahlberg and Ito 
[32], MATLAB® built-in function bvp4c implements a collocation method for the 
solution of the boundary value problem [33]. 

1.3.2 Constraints on the Low-Thrust Optimization Problem 

While solving the optimization problem, it is necessary to take in account the feasible 
limits of the low-thrust propulsion, i.e., to consider the maximum continuous thrust 
available by the propulsion method. Hence, it is required to apply the constraints 
on the optimization problem, which becomes effective once the optimal trajectory 
reaches the limits of the available specific thrust limit. 

Considering the optimal control problem stated eqn. (1.1) but now the control input 
u has the following constraints, 


\u\ < u 


_ max 


and the performance index is same as in eqn. (1.2), 


1 r*f 

J — - / u r Rudx 

2 jto 


(1.9) 


( 1 . 10 ) 



Then, according to the Pontryagin’s minimum principle, the optimal condition is 
given as, 


n(x*(t), \*{t),u*{t)) < 'H(x*(t),X*(t),u(t)) (1.11) 

It states [16] that of all control histories, u(t), corresponding to the neighborhood of 
the optimal trajectory, [x*(t), A*(f)], the one that minimizes the Hamiltonian is the 
optimal history, u*(t). 

1.3.3 Applicability of the Approximate Equations of Motion 
for Unconstrained Optimal Solution 

This thesis aims at examining the limits of the approximate equations of orbital 
motion, as stated in [3], so that, unconstrained optimization can be applied to it. For 
the purpose, different cases have been studied by varying the various parameters in 
the optimal control problem, such as, the time of the operation, the initial values of 
the relative distance and the relative velocity and the simulations have been carried 
out using the MATLAB® software. To check the accuracy of the numerical solution, 
it has been compared with the closed-form form solution obtained by analytically 
solving the approximate equations. Further, the effects of the oblateness of the 
planet are taken into account and several cases have been studied to understand the 
deviations in the optimal solution caused by the introduction of oblateness effect. 


1.4 Organization of the Thesis 

This thesis presents the study of the optimization of the low-thrust orbital missions. 
The organization of the thesis is as follows. Chapter 2 presents the review of the 
literature, mentioning the previous related work in the area. Chapter 3 presents the 
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mathematical model of the problem, starting with the simple case of the perfectly 
spherical Earth model, and then linearized the system to get the approximate HCW 
model. To make the model more realistic, effects of the oblateness of the planet 
Earth are also added in the model subsequently. Chapter 4 provides the optimization 
process for the different models discussed in the Chapter 3, including the closed-form 
solution for the approximated model. Chapter 5 presents the numerical analyses of 
the various optimization models discussed in the earlier chapters. It compares the 
results of the approximate and the non-linear equations of the relative motion for 
examining the safe applicability of the approximated model through various case 
studies of spacecraft rendezvous and stationkeeping missions. It further compares 
the errors which gets in the system when the effects of Earth’s oblateness are taken 
into account. Chapter 6 gives a summary of the results obtained in the research and 
finally the appendices contains the derivation steps involved in the various chapters, 
which were skipped for brevity. 



Chapter 2 


Review of the Literature 


Relative orbital transfer using low-thrust is a popular field of study due to its various 
applications and further possibilities in the field. Based on the type of the propul¬ 
sion system, the trajectory optimization can be classified as impulsive control, and 
continuous low-thrust control. Carter and Humi [28] investigates the fuel optimal 
rendezvous near a point in the general Keplerian orbit assuming bounded thrust and 
constant exhaust velocity. Schaub and Alfriend [25] and Vaddi et al. [26] discusses the 
different impulsive control methods for formation establishment and reconfiguration 
of satellites in the orbit. 

The sources of low-thrust are usually the unconventional propulsion techniques such 
as, the electric propulsion and the solar-sail systems. Since, the magnitude of the 
thrust obtained from such propulsion systems is relatively small [2], there is utmost 
requirement of an optimal trajectory for the maneuvering spacecraft. For the case of 
the low-thrust rendezvous, the objective of optimal control problems is to guide the 
maneuvering spacecraft from an initial position and velocity to a final state such that, 
the final position and velocity of the maneuvering spacecraft relative to the target is 
zero, while minimizing the fuel consumption, whereas, in the case of stationkeeping 
problem, there can be desired separation between the maneuvering and the target- 
spacecraft. On the basis of whether there are constraints on the applied input, since, 
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due to obvious reasons, the thrust generated by the low-thrust propulsion has limits 
on its magnitude, the optimization can be constrained or unconstrained. 

Euler [29] studied the optimal low thrust coplanar rendezvous between the maneu¬ 
vering and the target spacecraft using the approximate equations of motion. He 
compared the optimal trajectory solution obtained from the approximate equations 
with the optimal trajectory solution obtained using the non-linear exact equations 
obtained using the Newton-Raphson method. 

Hinz [22] studied the optimal low-thrust transfers between neighboring coplanar 
circular orbits using the approximate model [3]. The Clohessey-Wilshire model [3] 
approximates the equation of motion assuming circular reference orbit. Tschauner 
and Hernpel [5] presents the approximate model of the equations of relative mo¬ 
tion, however, elliptical reference orbits can also be accounted for. Marinescu [4] 
presents the minimum propellant optimal rendezvous of the two spacecraft without 
accounting for any constraints on the control input and considering the circular and 
the elliptical reference orbits, sequentially. For the optimization for the circular 
reference orbit, he made use of the approximate linearized equation of motion similar 
to Clohesey-Wiltshire model [3]. It is demonstrated in [4] that the maximum absolute 
values of the required control input falls in the domain of the low-thrust propulsion 
[2], Scott and Spencer [23] studies the optimal reconfguration of an n-spacecraft 
formation, considering the optimal low-thrust continuous transfer and making use 
of the Clohessy-Wiltshire model [3]. In their analyses, they have considered that 
the thrust input is unbounded, and given a set of transfer conditions, the optimal 
starting position for a single maneuver is derived. Based on the analytical solution 
to the performance index of a single transfer, this theory is extended to the full 
reconfiguration problem, which calls for a threefold optimization in terms of the 
dynamics and control of individual satellites, the order of reassignment, and the 
relative orientation of the final formation. 

Using the Tschauner and Hernpel model [5], Sengupta and Vadali [30] worked out 
the analytical solution for power-limited optimal minimum fuel rendezvous near an 
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elliptic orbit. Using the same model, Sharma et al. [31] presents an approach to the 
design of near-optimal feedback control laws, for minimum-fuel rendezvous between 
satellites in elliptic orbits of arbitrary eccentricity. The rendezvous problem for 
the nonlinear differential gravity model is solved by the application of neighboring 
optimal feedback control methodology used in conjunction with a nominal trajectory, 
obtained by solving the related minimum-fuel feedback control problem for the linear 
Tschauner-Hempel equations [5], analytically. 

Guelman and Aleshin [17] presents the bounded low-thrust, fixed-time, fuel-optimal 
constrained terminal approach direction rendezvous problem using the linearized 
equations of motion [3], in which they have developed a minimum-fuel two stage 
solution. The first stage of the solution, optimal transfer from the initial conditions 
to an intermediate point in the final line approach was implemented and in the 
second stage, an optimal thrust program guaranteeing directional approach was 
applied. Optimal rendezvous with fixed terminal approach is achieved by minimizing 
overall fuel consumption as a function of the intermediate point relative location and 
velocity in the final line approach. 

Since, due to approximation, the applications of the linearized models [3] and [5] 
are limited to cases when the separation between the maneuvering and the target- 
spacecraft- is relatively smaller as compared to the radius of the reference orbit-, 
non-linear model can be used to find the optimal trajectory. Shao et al. [24] studied 
the optimal fuel control design of the nonlinear spacecraft rendezvous system with 
the collision avoidance constraint, and formulated an optimal control problem with 
state-constraints by using exact penalty function method to transform the constrained 
opt-mizat-ion problem into a sequence of approximate unconstrained optimization 
problems. It was shown that the solution of these unconstrained problems converge 
to the solution of the original problem. 

Epcnoy [27] designed a numerical method for computing constrained fuel-optimal 
continuous-thrust rendezvous trajectories. By adapting the idea of using smoothed 
exact penalty functions to the optimal control context, it is shown that it is possible to 
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solve this kind of state-constrained problem through Potryagin’s minimum principle 
by simply solving a sequence of unconstrained problems. 

Wu et al. [18] developed the fuel-optimal low-thrust trajectories for high-precision 
formation maneuver using an exact nonlinear relative satellite dynamics model with 
orbital eccentricity and Earth oblateness effects. They also included collision avoid¬ 
ance constraints to ensure the safe maneuvering of the spacecraft. The resulting 
nonlinear optimal control problem is converted into a nonlinear programming problem 
by the application of the Legendre pscudospect.ral method and is then solved by 
using a sparse nonlinear optimization software. 

Till now, the major focus of research in the field of optimal low thrust relative orbital 
transfer is on the accurate modelling and formulation of the optimal control problem 
and obtaining optimal control using the approximate model. However, the applicabil¬ 
ity of the approximated model by varying the initial conditions of the maneuvering 
spacecraft, such as, the initial relative velcit.y and the initial relative distance, has 
not been examined. This thesis presents the first exhaustive numerical study of the 
optimal control problem of relative motion in the orbit using low-thrust. Different 
cases have been studied to account for the variations in the optimal control solution 
by varying the initial relative velocities of the maneuvering spacecraft, the initial 
separation between the maneuvering spacecraft and the target, and the specified 
time of the relative transfer. 



Chapter 3 


Orbital Relative Motion Dynamics 


In this chapter, the mathematical model for the relative motion in the orbit has 
been formulated, starting from the ideal case of a perfectly spherical planet, such 
that, there are no perturbation forces acting on the spacecraft, to a more realistic 
case, considering the oblat-cness of the planet Earth, and thereby introducing the 
perturbation forces into the model. 


3.1 Frames of Reference 

For analyzing the different relative motion models, introduction of the various frames 
of reference used is essential to be presented here. The various frames of reference 
used in the report are Earth-centered inertial frame (X), Earth-centered perifocal 
frame (V), local-vertical, local horizontal frame (£) centered at the spacecraft and 
an Earth-centered polar rotating frame {IV). The Earth-centered inertial frame is 
centered at the Earth, and the fundamental plane is the equator, x is directed from 
the Earth’s center to the vernal equinox, 5 is normal to the fundamental plane, 
positive towards the geographic north, and y completes the setup. The perifocal 
frame is centered at the Earth, and the fundamental plane is the orbital plane, x is 
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directed from the Earth’s center to the periapsis, z is nominal to the fundamental 
plane and is positive in the direction of the orbital angular momentum vector, and y 
completes the setup. Hill’s frame is a rotating frame of reference with respect to the 
Earth centered inertial frame and is centered at the chief. The fundamental plane 
for the Hill’s plane is the orbital plane, with x-axis directed from the spacecraft to 
radially outward, z-axis normal to the fundamental plane and positive in the direction 
of angular momentum vector, and y-axis completes the setup. The Earth-centered 
polar reference frame is centered at the Earth, and the fundamental plane is the 
orbital plane, r is directed radially outward from the Earth’s center and the angle 9 
is measured in the counterclockwise diretion from a refernce line. 

For the sake of brevity, throughout this work, the ordered-set notation is used 
to represent vectors, which shows the values of the respective components inside 
triangle brackets (( and )), separated by commas. For example, the position vector 
f — 2i + 9 j + 7k is represented as r = (2, 9, 7) in the ordered-set notation. 
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3.2 Equations of Relative Motion in the Orbit 

This section describes the non-linear equations of relative motion between the target 
spacecraft and the maneuvering spacecraft, which can use low-thrust to rendezvous 
or assume station-keeping with the target spacecraft. The spacecraft whose motion 
is independent of the other satellite, is known as the target, and is assumed to be in 
the circular orbit. The spacecraft which is chasing the target, and hence, needs to 
maneuver in a set pattern so as to maintain a desired distance from the target, is 
known as the maneuvering spacecraft. 

Consider two spacecraft in the orbits around the Earth, one of them being the target 
spacecraft which is in a circular orbit and the other one being the maneuvering 
spacecraft. Let f t denote the position vector of target spacecraft orbit, r m denote 
the position vector of maneuvering spacecraft orbit and yU = 398600.4 Km 3 /s 2 being 
the gravitational parameter of the Earth. 

For the target spacecraft, the equation of motion is, 

n + = 0 (3.1) 

' t 

and, similarly, for the maneuvering spacecraft, 

r m + ^~r t = 0 (3.2) 

^ m 

Now, if f denotes the position vector of the maneuvering spacecraft relative to 
the target spacecraft, directed from the target to the maneuvering spacecraft, it is 
obtained by subtracting the position vector of the target spacecraft from the position 
vector of the maneuvering spacecraft as follows: 


r = r m - r t 


( 3 . 3 ) 
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Figure 3.1: Configuration of Target and Maneuvering Spacecraft in the orbit: The 
position vectors of the target and maneuvering spacecraft are ft and f m , respectively. The 
position vector of the maneuvering spacecraft relative to the target spacecraft is f. The 
orbital frequency of the target orbit is n. The Hill’s frame is fixed on the target spacecraft. 
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Subtracting Eq. (3.2) from Eq. (3.1), we get, 


r t + r m 


= 0 


or, 


- b - 
r t = -~r t 


/x 


or, 


r — n 


r t 


I'm 

~3~ 


r. 


m 


from Eq. (3.3), r m = r + f t 


r — n 


t t 


r + r t 
|r + r t || 3 


(3.4) 


Eq. (3.4) is the relative acceleration between the target and the maneuvering space¬ 
craft. 


The angular velocity of Hill’s frame relative to the ECI frame is: 


0 = (o, 0. <j>) 


(3.5) 


and, position vector of the target, 


f t = (r t , 0, 0) 


(3.6) 


Let the coordinates of r in Hill’s frame be, 


f= (x, y, z) 


( 3 . 7 ) 
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Using Coriolis theorem, relative acceleration between the target and the maneuvering 
spacecraft in Hill’s frame is obtained as follows: 


- d 2 r _ dr du _ . _ . 

r = —- + 2u x — + — x r + u> x (u x r) 
dt z dt dt 

(3.8) 

Using the orbital reference frame, it follows that, 


(J 2 r 

dfi = <x - y ’ z) 

(3.9) 

dv 

dj X -^ = {0, 0, (j)) x (x, y, z) = (—(fry, <j>x, 0) 

(3.10) 

o' 

o' 

*3 14$ 

(3.11) 

^xr = (0, 0, 4>) x (x, y, z ) = (-(fry, 4>x , 0) 

(3.12) 

u x f = (—(fry, (frx, 0) 

(3.13) 

Lj x (u x r) = (—(fr 2 x, —<fr 2 y, 0) 

(3.14) 

Substituting Eqs. (3.9) - (3.14) into Eq. (3.8), and then separating the respective 

components along x-, y- and z-axis, we obtain, 

r x — x — 2 (fry — (fry — <fr 2 x 

(3.15) 

r y — y + 2(frx + (frx — (fr 2 y 

(3.16) 

r z = z 

(3.17) 


Resolving the respective components in Eq. (3.4), we get, 

1 x + r t 

x ^ r 2 [(x + r t ) 2 + y 2 + 2 2 ] 3 / 2 

_ y _ 

(x + r t ) 2 + y 2 + £ 2 ] 3 / 2 


(3.19) 



(3.18) 
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r z = y I - 


(x + r t ) 2 + y 2 + ^ 2 ] 3 / 2 


(3.20) 


From Eqs. (3.15)-(3.20), the non-linear equations of relative motion in the orbit are 
obtained as follows: 


x = 2 (f>y + cjyy + 0 2 x + y — 


x + r t 


r 2 [(x + r t ) 2 + y 2 + z 2 ] 3 / 2 


(3.21) 


y = —2 <fix — (j)x + <p 2 y + y 


(x + r t ) 2 + y 2 + z 2 ] 3 / 2 


z — y 


[(x + r t ) 2 + y 2 + ^ 2 ] 3 / 2 


(3.22) 


(3.23) 


Equations (3.21) - (3.23) constitute equations of relative motion in the orbit. 
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3.3 The Hill-Clohessy-Wiltshire Equations 

The HCW equations [3] presents a simplified and linearized model of relative motion 
of two spacecraft. 

For a circular orbit of the target spacecraft, we have a constant orbital frequency, 


d(f) 

— — n 
dt 

(3.24) 

where n is a constant. It follows that, 


0 = 0 

(3.25) 

also, 


r t = const. 

(3.26) 

Therefore, Eqs. (3.15) - (3.17) can be expressed as follows: 


f x = x — 2 ny — n 2 x 

(3.27) 

f y — y + 2 nx — n 2 y 

(3.28) 

■r z = z 

(3.29) 


Eqs. (3.18) - (3.20) result in the following: 

fy = A i 

f z = I 1 
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and, as n = \J / u/a 3 , Eq. (3.30) can be written as, 


2 3 


x + r t 


r t [(^ + r t ) 2 + y 2 + z 2 ] 3 / 2 


Using the approximation that relative distance between the target and the maneu¬ 
vering spacecraft is much smaller as compared with the radius of the reference orbit, 
y 2 + z 2 « (x + a) 2 , 

r i i 

(3.33) 


2 9 3 

r x = n r t — n a 


(X + r t ) 2 

and, by using Binomial theorem, Eq. (3.33) can be further simplified to the following: 


2 2 3 

r x = n r t — n r t 


1 2x 

rf n 


= n 2 r t — n 2 rt + 2xn 2 = 2xn 2 


(3.34) 


Similarly, we have, 


2 3 

r y = n r'l 


-y 


= -n y 


r z = —n z 


(3.35) 

(3.36) 


Comparing Eqs. (3.27) - (3.29) with Eqs. (3.34) - (1.36), we obtain: 


x = 3 nx + 2 ny 


(3.37) 


i) = —2 nx (3.38) 

z = — n 2 z (3.39) 

Equations (3.37) - (3.39) constitutes the linearized equations of relative motion for 
the maneuvering spacecraft. If there is an external control force acting per unit mass, 
u = (u x , u y , u z ), we have above equations modified as: 


x = 3 n 2 x + 2 ny + u x 


(3.40) 


y = —2 nx + u y 


( 3 . 41 ) 
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z = — n 2 z + u z (3.42) 

However, the model formulated by Clohessy and Wiltshire [3], is very simple, but 
is limited to application by the assumptions made in the process of linearization, 
such as the results starts varying when the relative separation between the target 
and the maneuvering spacecraft increases beyond a certain safe applicability limit, 
which is also discussed in the present work. Furthermore, the model assumes that 
the target spacecraft is in a circular orbit around a spherical planet, rendering it 
useless if the orbit of the target is not circular, for example if it is elliptical. Since, 
at comparatively low-Earth orbits, there is considerable resistance offered by the 
atmosphere, the above model does not give accurate results. 

To account for the eccentric orbits of the target, Tschauner and Hempel [5] developed 
a model for the relative motion of the spacecrafts in the orbit, but instead of treating 
time as the independent variable, as in the case of Clohessy and Wiltshire [3], they 
took true anomaly of the reference target orbit as the independent variable, the result 
being three, second-order linear differential equations. Unlike the model developed by 
Clohessy and Wiltshire [3], the Tshauner and Hempel model [5] for relative motion 
in the orbit can also take problems in which the target is in an elliptical orbit. 
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3.4 State Space Representation of the HCW Model 


For further analyses, it is necessary to convert the above model to its state-space 
form. 

Let us define the state variables to be, 

x\ = x, X 2 = x, x 3 = y, £4 = y, x 3 = z and x§ — z 

such that, x = (xi, x 2 , x 3 , x 4 , x 5 , x 6 ) 

and, 

U (u XI U yi U z ) 

Hence, the HCW model derived in section (3.3) can be written in the state-space 
form as x = Ax + Bu 


X\ 


0 

1 

0 

0 

0 

0 


X\ 


0 

0 

0 

X 2 


3 n 2 

0 

0 

2 n 

0 

0 


x 2 


1 

0 

0 

X3 

_ 

0 

0 

0 

1 

0 

0 


X3 

+ 

0 

0 

0 

X\ 


0 

—2 n 

0 

0 

0 

0 




0 

1 

0 

x 5 


0 

0 

0 

0 

0 

1 


x 5 


0 

0 

0 

x 6 


0 

0 

0 

0 

— n 2 

0 


x 6 


0 

0 

1 


^ X 

u y 

u z 


(3.43) 


The linearized state-space formulation will be used for carrying-out the optimal 
control calculations in the next chapter. 
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3.5 Effects of Orbital Perturbations 

To make the relative motion model more realistic and useful for practical application, 
it is necessary to consider the perturbations in the orbit caused by various factors, such 
as due to atmospheric drag and the oblateness of the planet. Effects of atmospheric 
drag are negligible in high-Earth orbits, but if we are considering low-thrust maneuvers 
in the low-Earth orbits, the atmospheric drags becomes considerable. Carter and 
Humi [6] formulated the orbital relative motion equations taken in account the effects 
of quadratic drag and with an assumption that the reference target orbit is not highly 
eccentric. 

Inclusion of the effects of the zonal harmonics of the Earth makes the model a little 
bit complicated, since, the orbital elements, such as inclination, argument of perigee 
and ascension of the node changes with time in the presence of orbital perturbations. 
The perturbation effects due to the oblateness of the planet Earth, J 2 , which is also 
discussed in the present work, are quite significant and causes the precession of the 
orbital frame. Kechichian [7] formulated the exact nonlinear model of relative motion 
in the orbit taking in account the effects of perturbation due to J 2 and atmospheric 
drag. The model, being very complex and non-linear, makes it difficult to employ in 
control applications. 

Schweighart and Sedwick [8] presents a high-fidelity linearized J 2 model for the 
satellite formation flight, taking the time-average of the gradient of the J 2 potential 
to calculate the in-plane relative motion of the spacecraft for a circular reference 
orbit. There are several advantages linked to this model, since, it consists of linear, 
constant-coefficient differential equations, which can easily provide analytic solutions. 
Furthermore, the model also discusses the effects of the period mismatch between 
the in-track and cross-track motion, which then produces a tumbling effect and the 
differential drift which causes the break-up of the satellite formation over time. 
Sabatini et al. [12] discusses the modelling of relative motion considering the 
perturbed environment. They have used an alternative approach, which has also 
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been used by Izzo and Valente [13], to include the J 2 differential effects in a linear 
model. Starting with the reference Keplerian orbit, they obtained the expression for 
the radius of the perturbed reference orbit, which in-turn is used for deriving the 
linearized equations of relative motion in the orbit. 


3.6 Model Relative to Perturbed Reference Orbit 

In order to derive the equations of motion relative to the perturbed reference orbit, 
it is required to find the expression for the radius for the perturbed reference (target) 
orbit. For this purpose, it is required to evaluate the difference in the target’s orbit if it 
is purely Keplerian from when it is having the perturbation due to Earth’s oblateness, 
similar to the approach used by Izzo and Valente [13]. Let the perturbation forces 
be represented by F p , then, the two-body equations for the purely Keplerian orbit 
and the one with perturbation forces present can be expressed as follows: 

r k = -’4 r ~lc (3.44) 

n 

' E p = + F p (3.45) 

p 

where, the subscripts k and p represents purely Keplerian and the perturbed orbit, 
respectively. Now, the perturbation force, F p is given as (Appendix A), 

F P = ( f r , fe, fh > 

3 J uH? 

= - 2 , 6 (1 — 3 sin 2 i sin 2 6, sin 2 i sin 2d, sin 2 i sin 0) (3.46) 

2 Tq 

For deriving the changes in the radius of the target spacecraft due to planatery 
oblateness, it is required to take the difference of the Eqs. (3.45) and (3.44): 

r e = r p - r k 


( 3 . 47 ) 
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where, r e = (x e , y e , z e ) 

Thus, the following set of differential equations are derived: 


x e = 2 ny + 3 n 2 x e + rj(l — 3 sin 2 i sin 2 9) 


(3.48) 


ij e = —2 nx e + r/(sin 2 i sin 29) 

(3.49) 

z e = —n 2 z e + rj (sin 2 i sin 9) 

(3.50) 


where, 


V = 


3 J 2 yU\ 


2 r 0 4 


(3.51) 


The linearized differential equations Eqs. (3.48)-(3.50) have the particular solution 


for XX. 


Xe 


3 J 2 R\ 
2 To 




(3.52) 


Thus, the expression for the radius of the reference orbit is obtained by adding ro to 
the particular solution obtained in Eq. (3.52): 


r = r 0 + x e 

, 3 J 2 K\ 

(3.53) with respect to time, we derive: 

3 J 2 nR 2 e 
2 r o 

The orbital angular momentum is defined as follows: 



2 r 0 

Differentiating Eq. 


r = 


- sin" i cos" 6 H— sin" i — 1 + ( 1-sin" i cos 9 


(3.53) 


(3.54) 


h = r x v = r x (r + lJ x r) 


( 3 . 55 ) 



As the orbital plane is defined by the position and velocity, the velocity of the target 
spacecraft must be in the reference frame, so that, f is parallel to r, 

h — rx{Cjxr) (3.56) 


Using the property of the vector product 1 , Eq. (3.56) can be expressed as: 


h — u ■ r 2 — r(r ■ to) 


(3.57) 


Rearranging, 


^. h ? 
u — (r • u)r H —-r n 


(3.58) 


Evaluating the time-derivative of angular momentum, 


dh d 
— — — (r x v) 
dt dt K 

— fxv + fxv 


= v x V + r x (-- + Fp) 


= f x F n 


(3.59) 


where, F p can be expressed in its radial, tangential and normal components as in 
equation (3.46), 

F p = (. f r , fe, U) (3.60) 

so that, 

dh 

- = <r, 0, 0) x (f r , f e , f h ) 

= ( 0 , rf h , rf 9 ) 

= rf h 9 + r f g h (3.61) 


1 Property of Cross-Product: a x (b x c) = b(a ■ c) — c(a ■ b) 
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The above time-derivative of angular momentum can also be expressed as, 


dh ? _ r 

— = h + to x h 
at 


(3.62) 


Substituting the value of uj from Eq. (3.58) into Eq. (3.62), the following expression 
is obtained for the time rate of change of angular momentum: 


dh b {h ^ , 

~T 7 = h + \ —; + r(r ■ uj) ] x h 


dt \r" 

= h + r (r • c u) x h 


(3.63) 


Since, it is true that, 


6 — r x h 


(3.64) 


Eq. (3.63) can be written as follows: 


^ = /?/?. + (f • u)hO 


(3.65) 


Comparing Eqs. (3.61) and (3.65), we obtain: 


h = rf e 


(3.66) 


(f -uj)h = rf h 


to get, 


r ■ u = 


rfh 

h 


(3.67) 


Substituting Eqs. (3.46) and (3.67) into Eq. (3.58), we obtain the following result 
for the orbital angular velocity: 


u = 


3 J2trR 2 e 
2 r 3 h 


sin2isin#, 0, 



(3.68) 
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By differentiating Eq. (3.68) with respect to time, we obtain: 

, 0, 


3 T ^. / nrh cos 9 (rh + 3hr) sin 9 \ „ / hr — 2hr 


co = (—-J 2 nR e sin 2i 

2 \ 


r A h 2 


I) (3.69) 


The orbital motion of a spacecraft orbiting around a planet is given by the two-body 
equation as follows: 

it 2 r n 

(3.70) 


d 2 r a _ 

7.o o 1 I ± V 


dt 2 

For the target and the maneuvering spacecraft, with the oblateness effect taken in 
account, the respective equations of motion are as follows: 


d 2 f t 
dt 2 

rpr 

LL I m 

dt 2 


3 




h + F pt 

(3.71) 

Tm + F pm 

(3.72) 


By Taylor’s series expansion of Eq. (3.72) around the target’s orbit, the following 
result is obtained: 


d 2 f 

U, I m 

dt 2 


H „ / H . 

~^r t + V r t 

rf W 


+ V —»r t 

n \n 


(■ r m - r t ) 


+ Fp t + V F pt ■ (r m — rt) 


+ I Fpt + V.Fpt • f 


(3.73) 


Subtracting Eq. (3.70) from Eq. (3.73), and proceeding in the same manner as in the 
unperturbed case, we get the position vector of the maneuvering spacecraft relative 
to the target spacecraft, 


+ 2cJx y + f + to x (cJ x r) = —V f ^f t ) • r + WF P ■ r (3.74) 

dt 2 dt dt \ r t 


In the above equation we have two vector functions given in the spherical coordinates, 
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viz., 


gf) 


0 ,0> 


(3.75) 


and 


F pt = — ^ (1 — 3 sin 2 i sin 2 9, sin 2 i sin 29, sin 2i sin 9) 


2 rf 


(3.76) 


The respective tensors [ 8 ] are calculated as (Appendix B), 


V g(r) = 


2f 


0 

0 


JL 

r 3 


0 

0 


(3.77) 


and, 


VF Pt = 


6J 2 yR i 


1 — 3 sin 2 i sin 2 9 
sin 2 i sin 29 

sin 2% sin 9 


sin 2 i sin 29 

■\ + sin 2 if l sin 2 9 - \ 


— | sin 2 i cos 9 


sin 2 i sin 9 


-1 sin 2 i cos 9 


-§ + sin 2 i I | sin 2 0 + | 


(3.78) 


Substituting the values of tensors from equations (3.77) and (3.78) in equation (3.74) 
and expanding it in terms of the respective components in the Hill’s frame, the 
following equations are obtained: 


2 p 6,/2/iA 2 22 

x = 2u z y + oj z x — uj x u z z + (j o z y + 2—^x H- g —(1 — 3 sm 1 sin 9)x 


rf 


rf 


6J 2 nR 2 el . 2 , . om „. , 6 J 2 yR 2 e 


(sin 2 i sin 2 9)y + 


/v> 5 
' t 


(sin 2i sin 9)z 


- x\ to z + 2-3 + 

' t 


H 6J2I-1R') 


(1 — 3 sin 2 i sin 2 9) ) + y ( cu z + 


6J 2 yR i 


•2 

e /„■ 2 


sin 2 i sin 29) 


6J 2 yRl . . 

+ z [ — u) x u> z H-g— (sm 2i sm 9) ) + y(2u z ) 


(3.79) 
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2,2 • I • A 4 I 6,/ 2 /J -R e / ■ 2 


y = — 2 u z x + 2u; x i + t o z y + out/ - cu z a: + d; x z-o y + 


sin 2 isin2d)a; 


QJ 2 fiR 2 e 


7 . 1 

-h sin" 1 i I - sin 2 0- 


J/ + 


6 


e 1 — - sin 2i cos 9 ) z 


= x 


. 6 J 2 /ii? 2 . 2 . . 

H-=■— sm i sin It) 


2,2 F I 6J 2 /i/?, 
+ W * “ 73 + 15 


1 . 2 .,7 . 2/1 1, 

- +sm i(-sm 0 - - 




QJ 2 yR 2 e f l sin2icos0 j +u 2 


(3.80) 


.. 2 . u 6 J 2 uR 

z = -2 u x y - uj x uj z x + uj x z - uj x y - z + 


ri 


n 


QJ 2 jiR 2 e 


— - sin2 i cos 6 )y + 


= x 


+ z 


w x w z + 


6J 2 /ii? 2 


77 


6 J 2 fiRl 


(sin 2i sin 0) ) + y 


sin2i sin0)x 


5 


2 H . b-hi-iR 2 

(jJ~ — o I 

/V >0 

•t 


<y>5 

' t 


4 

+ 


3 • 2-/5 . 2 n 1 

- ^ +sm *(-sin 0 + - 


—- + sin i - sin 9 + - 


4 

6 J 2 yR 


- sin 2* cos 0 
4 


4“ i/( 2 u x ') 


(3.81) 


For further analyses and for solving the optimal control problem, the above model 
is necessary to be represented in the state-space form. Same as in the case of the 
HCW model, the state-variables are defined as, 
x\ = x, x 2 = x, X 3 = y, X 4 = ij, X 5 = z and x§ = i 

Hence, the model derived in Eqs. (3.79) - (3.81) can be represented in the state-space 
form as, 


x\ 


0 

1 

0 

0 

0 

0 


Xl 


0 
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0 

x 2 


9i 

0 

92 

2u z 

93 

0 


X 2 


1 

0 

0 

X3 

— 

0 

0 

0 

1 

0 

0 


X3 

+ 

0 

0 

0 

X4 


94 

-2 u z 

95 

0 

96 

2(jj x 


X4 


0 

1 

0 

x 5 


0 

0 

0 

0 

0 

1 


x 5 


0 

0 

0 

x 6 


97 

0 

98 

—2 c o x 

99 

0 


x 6 


0 

0 

1 


u y 

u z 


(3.82) 
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where, 


and, 


qi = uj 2 z + 2n 2 + K (1 — 3 sin 2 i sin 2 9) 
q 2 — u z + A' (sin 2 i sin 29) 
q 3 = —(jj x u) z + K (sin 2 i sin 9) 
g 4 = — oj z + K (sin 2 i sin 29) 


2 i 2 2 i tv" 

+ U x ~ n + K 


l 


7 


-h sin i - sin 9 - 


g 6 = 6 j x — -A' sin 2i sin d 
qi = — oj x oj z + K sin 2i sin d 


q 9 = col - n 2 + A' 


g 8 = — ui x — -K sin 2 i cos 9 

3 . 2 ./5 . 2/1 1 

-h sin i - sin d H— 


(3.83) 

(3.84) 

(3.85) 

(3.86) 

(3.87) 

(3.88) 

(3.89) 

(3.90) 

(3.91) 


C././( /i 1 ) 

<y>5 
' t 


(3.92) 


Eq. (3.82) gives the state-space model for HCW equations including the effects of 
Earth’s oblateness effect. 


3.7 Modelling of Relative Motion Using Orbital 
Elements 

Hamel and Lafontaine [9] presents the linearized dynamical model of formation flying 
spacecraft on a J 2 perturbed elliptical orbit. They have developed an analytical 
state-transition matrix that models the relative motion about the elliptic reference 
orbits under J 2 perturbation, so that, it can be solved to give analytical results. 
They have used the geometric approach similar to as used by Gim and Alfriend [10], 
but their model neglects the variation in the short-period relative motion induced by 
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the J 2 perturbation between the target and the maneuvering spacecraft, however, 
includes the relative mean orbit element- drift to the natural osculating element 
Keplerian dynamics. 

Alfriend et al. [11] discusses the modelling of the relative satellite motion using orbital 
elements. The main advantage of using the orbital elements in the equations of the 
relative motion is that they allows the straightforward incorporation of the orbital 
perturbations, yielding Lagrange’s planetary equations and the Gauss variational 
equations. It also facilitates the derivation of higher order nonlinear extensions to 
the HCW solution in the Cartesian coordinates. Furthermore, it also models the 
perturbed relative motion using orbital elements, and thereby, deriving the linearized 
^-differential equations for circular orbits. It discusses the formation control using 
discrete-time LQR method and the simulation of the formation flying. 



Chapter 4 


Optimization of the Low-Thrust 
Motion 


The following chapter illustrates the optimization process. The objective here to 
design a fuel-optimal low-thrust control input to be applied to the system for the 
desired operation in the orbit. The optimal control approach has been adapted 
from [16]. The various operations considered here are coplanar rendezvous of the 
maneuvering spacecraft with the target, rendezvous of target with maneuvering 
spacecraft taking in account the out-of-plane motion also, and the same has been 
repeated with the effects of Earth’s oblateness taken in account. 
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4.1 Optimization of the Control Input Trajectory 
for a General Case including the Effects of 
Earth’s Oblateness 


Consider the maneuvering spacecraft is at some distance, r 0 = (xq, yo, z 0 ) and 
velocity, Vo = ( v Xo , v yo , v Zo } relative to the target spacecraft, which is in a circular 
orbit. The thrust can be applied along all the three axes as u = ( u x , u y , u z ). 
Considering the relative motion model which takes in account the effects of Earth’s 
oblateness developed in section 3.7, we have the state-space model obtained as 
follows: 


Xi 


0 

1 

0 

0 

0 

0 


X\ 


0 

0 

0 


X2 


Qi 

0 

92 

2ui z 

93 

0 




1 

0 

0 


r -| 
















^X 

x-s 

= 

0 

0 

0 

1 

0 

0 


X3 

+ 

0 

0 

0 


ILy 

x 4 


94 

-2c o z 

95 

0 

96 

2u x 


X 4 


0 

1 

0 

















u z 

X5 


0 

0 

0 

0 

0 

1 


X 5 


0 

0 

0 



X 6 


97 

0 

98 

—2 u x 

99 

0 


x 6 


0 

0 

1 


the parameters viz. 

to 

kQ 

CO 

94, 

95, 96: 

. 97 , 9s and q 9 

are defined 


(4.1) 


3.7. 


Here, the interest is to design a fuel optimal control system. For fuel-optimization, 
the above model is subjected to performance index, 


J 


1 

2 

1 

2 


'to 


( u T Ru)dt 


( u l + u l + u 2 z )dt 


to 


(4.2) 
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The Hamiltonian for the above system is obatined as follows: 

% = u T Ru + A T f 


(4.3) 


where, 


/ = 


0 

1 

0 

0 

0 

0 


aq 


0 

0 

0 

qi 

0 

92 

2 u z 

93 

0 


(7'2 


1 

0 

0 

0 

0 

0 

1 

0 

0 


x 3 

+ 

0 

0 

0 

94 

- 2 u z 

95 

0 

96 

2 u x 


X 4 


0 

1 

0 

0 

0 

0 

0 

0 

1 


x 5 


0 

0 

0 

97 

0 

98 

— 2 uj x 

99 

0 


X 6 


0 

0 

1 


^X 

lly 

u z 


(4.4) 


Therefore, Eq. (4.3) can be written as follows: 

111 

R = -u^. + -Uy + -U 2 Z + X\X 2 + X2{q\X\ + q 2 Xz + + Q3 x 5) + A3(:r 4 ) 

+A 4 (g 4 aq - 2u z x 2 + q 5 X 3 + q 6 x 5 + 2uq£ 6 ) + A 5 a; 6 + A 6 (g 7 a;i + q$x 3 - 2w x x 4 + g 9 a; 5 ) 

(4.5) 

As per the necessary condition for optimality, we obtain: 


A* = 


dR 

dx 


= —AX 


T \ * 


and here, 


A = 


q i 
o 

94 

o 

9? 


1 0 0 0 0 

0 q 2 2 u z q-z 0 

0 0 10 0 
-2u z q 5 0 g 6 2 u x 

0 0 0 0 1 


0 


98 


-2u x q 9 0 


(4.6) 
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The condition for stationarity (necessary condition for optimality) is as follows: 


dU 

du 


= 0 


which results in, 


* \ * 

= -A 2 

><;, = -a: 

< = -a; 


(4.7) 


(4.8) 

(4.9) 
(4.10) 


Clearly, the system forms a two-point, boundary value problem with time-variant 
terms in the state-matrix, ft. can be solved numerically to give the optimal input, 
and is discussed in Chapter 5. 


4.2 Optimization of the Control Input Trajectory 
without Oblateness Effects: Closed-Form So¬ 
lution 

Consider the maneuvering spacecraft is at some distance, r 0 = (xo, yo, Zo) and 
velocity, v 0 = ( v XQ , v yo , v ZQ ) relative to the target spacecraft, which is in a circular 
orbit. The thrust can be applied along all the three axes as u = ( u x , u y , u z ). 
Considering the case when there is no perturbation due to oblateness effects of the 
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planet, the linearized state-space model reduced to that mentioned in Eq. (3.43). 




0 

1 

0 

0 

0 

0 


X\ 


0 

0 

0 

^2 


3 n 2 

0 

0 

2 n 

0 

0 


x 2 


1 

0 

0 

^3 

_ 

0 

0 

0 

1 

0 

0 


X3 

+ 

0 

0 

0 

0:4 


0 

—2 n 

0 

0 

0 

0 


X4 


0 

1 

0 



0 

0 

0 

0 

0 

1 


x 5 


0 

0 

0 

X(i 


0 

0 

0 

0 

— n 2 

0 


x 6 


0 

0 

1 


^X 

u y 

u 7 


(4.11) 


subjected to performance index, 


1 f tf 

J = - / (u 1 Ru)dt 

2 Jt 0 

= \J t (“»+ “»+ u D dt 


(4.12) 


The Hamiltonian is obtained as follows: 


TL = u T Ru + A T f 


(4.13) 


where, 


/ = 


0 

1 

0 

0 

0 

0 


X\ 


0 

0 

0 

3 n 2 
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0 

2 n 

0 

0 


X 2 
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0 

0 

0 

0 

0 

1 

0 

0 


X 3 

+ 

0 

0 

0 

0 

—2 n 

0 

0 

0 

0 


£4 


0 

1 

0 

0 

0 

0 

0 

0 

1 


X 5 


0 

0 

0 

0 

0 

0 

0 

— n 2 

0 


x 6 


0 

0 

1 


Therefore, Eq. (4.13) can be written as: 


^X 

u y 

u 7 


(4.14) 


1 1 1 
1 .2 i ± „,2 i 1 2 


T-L = ~u x + -u v + -u + Aix 2 + A 2 (3n xi + 2 nr 4 + u x ) + A 3 (a; 4 ) 
2 2 y 2 


+^(— 2nX2 + w„) + A 5 X 6 + Ag(— n^x 5 + u z ) 


(4.15) 
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As per the necessary condition for optimality: 


and here, 


A* = - 


A = 


m 

dx 


= —A A 


T \ * 


o 

3 n 2 

0 

0 

0 

0 


1 0 0 

0 0 2 n 

0 0 1 

-2 n 0 0 

0 0 0 

0 0 0 


0 0 
0 0 
0 0 
0 0 
0 1 
— n 2 0 


(4.16) 


Assuming initial values for A (t) as: A T (0) = [Ai(0), A 2 (0), A 3 (0), A 4 (0), A 5 (0), A 6 (0)] = 
C = [Ci, C 2 , C 3 , C 4 , C 5 , C 6 ], and as, 


A*(t) =e~ AT, C 


(4.17) 


We have, 


/ 


„ -A T t = 

4 — 3 cos (nt) 

sir(n t) 
n 

0 

2 2 co t) 

n n 

0 

0 


3 n sin (nt ) 6n t — 6 sin(nt) 6n cos(n t) — 6n 
cos (nt) 2 co ^ n ^ — ^ 2 sin(nt) 

0 1 0 
3 1- 

0 
0 


—2 sin(nt) 
0 
0 


4 sir(n t) 
n 


4 cos (nt) — 3 
0 
0 


0 

0 

0 

0 


0 

0 

0 

0 


\ 


cos (nt) n sin (nt) 
-liayl cce(ni) I 


(4.18) 
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Therefore, Eq. (4.17) results in: 


A *(t) = e~ AT C 

—Ci (3 cos (nt) — 4) — C4 (6 n — 6 n cos(nt)) — C3 (6 sin(nt) —6 nt) — 3 C2 n sin(nt) 
C 2 cos (nt) - C 3 (j - + 2 C 4 sin (nt) - Cl s ^ nt) 


C 3 

Cl (1 - 2 -^) + C 4 (4 cos (nt) - 3) - 2C 2 sin (nt) + C 3 ( 3 t- 

C 5 cos (nt) + C 6 n sin (nt) 

C e co, 

(4.19) 


And, for stationarity (necessary condition for optimality), 


which gives, 


dU 

du 


= 0 


u n 


u. 


u„ 


-a; 

-a: 


(4.20) 


(4.21) 

(4.22) 

(4.23) 
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Hence, we obtain: 


u* x = C, ( - 


2 2 cos(nt) 

n n 


C 2 cos(nt) — 2 C 4 sin(nf) + 


Ci sin(nf) 


n 


(4.24) 


u 


* 

y 


2 C 2 sin(nt) — C 4 (4 cos (nt) — 3) — Ci 

-c 3 



2 


4 


cos (nt) 


sin (nt)\ 

n J 


(4.25) 


u 


* 

z 


C 5 sin (nt) 
n 


Cg cos (nt) 


And as, 


where, 


x(t) = e At x(0) + / e A{t - T) Bu(r)dr 


B 


0 0 0 
1 0 0 
0 0 0 
0 1 0 
0 0 0 
0 0 1 


If we consider the initial states of the state-space model (4.11) as, 


(4.26) 


(4.27) 


(4.28) 


x(0) = [x 1 (0),x 2 (0),X3(0),X4(0),X5(0),x 6 (0)] t (4.29) 
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Now, 


/ 


e At = 


4 — 3 cos(nt) ^ 0 

3 n sin (nt) cos (nt) 0 

6 sin(nt) — 6 nt - 1 i iA4zA)_ 3t 

6 n cos (nt) — 6 n —2 sin (nt) 0 4 cos (nt) — 3 

0 0 0 0 

0 0 0 0 


2 2 cot) 

n n 

2 sin(nt) 


(4.30) 


From Eqs. (4.29) and (4.30), we obtain: 


0 

0 

0 

0 

cos (nt) 


0 

0 

0 

0 

sir(n t) 




-n sin (nt) cos (nt) J 


e At x(0 ) = 


/ 


*4(0) (I - 2 -^) - *4(0) (3 cos(nt) - 4) + 

x 2 (0) cos (nt) + 2x 4 (0) sin (nt) + 3nx 4 (0) sin (nt) 
x 3 (0) — x 2 (0) (f ~ 2 co ^ nt) ^j +xi(0 ) (6 sin(nt) - 6nt) - x 4 (0) (3t- 


\ 


4 sir(n t) 


x 4 (0) (4 cos (nt) — 3) — x 4 (0) (6n — 6n cos (nt)) — 2x 2 (0) sin (nt) 
x 5 (0) cos(nt) + Mo) ^ nt) 
x 6 (0) cos (nt) — nx 5 (0) sin (nt) 


(4.31) 


And, 
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pA(t—r) _ 


4-3 (9 


i e o J(i -1?) o o 


o o 


6g — 6n(t — r ) -(i9 •— 1 ) 1 3r — 3t + -g 0 0 


6 n($ — 1) 


-2 g 0 4(9-3 


where, i9 = cos (n (t — r)) and g = sin(n (t — r)). 
Multiplying Eqs. (4.32) and (4.28), we obtain: 


0 0 


(9 -g 

n ^ 

—ng (9 


(4.32) 


si r(n ( t—r )) 


2 2 co (t—r)) 


cos(n (t — r)) 


2 sin(n (t — r)) 


e A{t ~ T) B = 


2 co^n (t—r)) 2 g g ^ _i_ 4 sic(n (t—r)) 

n n n 

—2 sin(n (t — r)) 4 cos(n (t — r)) — 3 


sii^n (t—r)) 


cos(n (t — r)) 


Now, from Eq. (4.33) and noting that u = (u x ,u y ,u z ), we obtain: 


(4.33) 


Mq-r) fi ^(r) = 


66 + “66 
66 - 266 

66 — 6 [3t — 3t + ^6] 
-266-6(46-3) 
466 
-66 


(4.34) 
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where, 


£1 


— C\ (1 — cos (nr )) + C 4 (4 cos (nr) — 3) — 2C * 2 sin(nr) + C * 3 
n 


3 r-sin(nr) 

n 


£9 = —6*3(1 — cos (nr)) — C*9Cos(nr) — 2C*4sin(nr) H—Cisin(nr) 
n n 

£3 = sin (n(t - t)) 

1 

£4 = C* 6 cos (nr) - C* 5 sin(nr) 

n 

& = -(& - 1 ) 

n 

= cos(n(t - r)) 


(4.35) 

(4.36) 

(4.37) 

(4.38) 

(4.39) 

(4.40) 


Now, integrating Eq. (4.34) from 0 to t, we get, 


( 6,6 + -^ 3 ) dr 


4 C 2 16 C 3 
n 2 n 3 

11 C 4 sin (nt) 


+ 


4C 2 cos(nt) I 6 C 3 cos(nt) 13 Ci sm(nt) 
n 2 n 3 2 n 3 

3C 3 t 2 4 C] t 6 C 4 t 5Citcos(nt) 

n n 2 n 2 n 2 

5 C 4 t cos (n t) 5 C 2 1 sin(n t) 5 C 3 1 sin(n t) 

n 2 n n 2 


(4.41) 


f (tee ~2Us) dr 
Jo 

= -^-^(8 Ci — 12 C 4 n — 8 Ci cos (nt) — 22 C 3 sin (nt) + 12 C 3 n t 
+ 12 C 4 n cos (n t) — 3 C 2 n sin(n t) + 10 C 3 n t cos (n t) 

— 5 Ci nt sin(n t) + 5 C 2 n 2 1 cos (nt) + 10 C 4 n 2 1 sin(n t )) (4.42) 
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(66 - 6 [3 T -3 t + - 6 ]) dr 
i n 


28 C 4 


n z 


16 Ci 3 C 3 1 3 

-- H--— 

n 3 2 


9 C 4 t 2 16Cicos(nt) 28C 4 cos(nf) 

I n o 


llC2sin(nt) 38C3sin(nf) 

o I n I 


rr 


rr 


2 

3Cit 2 

n 


rr 


rr 


6 C 2 f 28 C 3 1 5C’2t cos(nt) 


n 


rr 


n 


10 C 3 1 cos(n t) 5 Ci t sin(n t) 

O I o 


10 C 4 t sin(nf) 


rr 


rr 


n 


(4.43) 


/'(- 266 - 6 ( 16 - 3 ))^ 

Jo 

9 C 3 1 2 6 C 2 28 C 3 6 C 2 Cos(nt) 28 C 3 cos (nf) llCisin(nf) 

: —n- 2 -9 C 4 1 H-1-n- 2 - 

2 n n n rr rr 

18C 4 sin (nt) __ . . __ . . , 6 Ci £ 

H- 10 C 4 1 cos (n t) + 5 C 2 1 sin(n t) 4 - 

n n 

5Cit cos(nt) 10 C 3 f sm(nt) ,, 
n n 



C 5 


f t co£) 


n 


2 


sir(n t) \ 
2 n y 


C 6 t sin(nt) 
2 n 


(4.45) 


and, 



66) dr 


Csf sin(nt) 
2 n 


C 6 


£ cos (nt) 


+ 


sin(nt) ^ 
2 n ) 


(4.46) 


The complete solution of the states is given by Eq.(4.27). Incorporating the solutions 
given by Eq.(4.31) and Eq.(4.41)-(4.46) in Eq.(4.27), the separate state-solutions are 
obtained as: 
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. . . . (2 2 cos (nt)\ , . , . x 2 ( 0 ) sin (nt) 

\n n J n 

4C 2 I 6 C 3 4C 2 cos (nt) I 6 C 3 COS (nt) 13Cisin(nt) 

4~ o 4” Q O O + 


n n u 
11 C 4 sin (nt) 


n 


rr 


2 n 3 


3 C 3 1 2 


4Ci t 6 C it 5Citcos(nt) 
O 4~ 


n 


n 


n 


n 


2 n 2 


5 C 4 1 cos (n t ) 5 C 2 1 sin(n t) 5 C 3 1 sin(n 1 ) 


n 


2 n 


n z 


(4.47) 


x 2 (i) = x 2 (0) cos (nt) + 2 x 4 (0) sin(nt) + 3nxi(0) sin(ni) 

— —^ (8 Ci — 12 C 4 n — 8 Ci cos(n 1) — 22 C 3 sin(n t) + 12 C 3 n t 
+ 12 C 4 n cos (n t) — 3 C 2 n sin(n t) + 10 C 3 n t cos (n t ) 

— 5 Ci nt sin (nt) + 5 C 2 n 2 1 cos (nt) + 10 C 4 n 2 1 sin (nt) ) (4.48) 


. . . . ... 2 2 cos (nt) 

(t) = x 3 ( 0 ) — x 2 ( 0 )- 


n 


n 


+ xi( 0 ) (6 sin (nt) — 6 nt) 


x 4 ( 0 ) ( 3 1 


4 sin (nt) 


n 


28 C 4 16 Ci 3 C 3 1 3 9 C 4 1 2 I 6 C 1 cos (nt) 

4 ~ o o 4 ~ 4 ” o 


rr 


rr 


rr 


28C 4 cos (nt) llC 2 sin(nt) 38C 3 sin(nt) 3Cil 2 

O I .1 I o + 


6 C 2 1 


rr 




rr 


n 


n 


28 C 2 ,t 5C 2 tcos(nt) 10C 3 tcos(nt) | 5Ci t sm(nt) 10C 4 tsin(ni) 

o o I o 


rr 


n 


n z 


n z 


n 


(4.49) 
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X 4 (t) — x 4 (0) (4 cos (nt) —3) — xi(0) ( 6 n — 6 n cos (nt)) — 2x 2 (0) sin (nt) 

9 C 3 1 2 6 C 2 28 C 3 6 C 2 cos (nt) 28 C 3 cos (nt) llCisin(nt) 

H---o-9 C 4 1 H-t-o- 2 - 

2 n n z n n z n z 

18C 4 sin (nt) .. __ . . , 6 Ci £ 

H-10 C 4 1 cos (n t) + 5 C 2 1 sinfn t) H- 

n n 

5Ci£cos(n£) 10 C 3 £ sm(nt) 

H-1-(4.50) 

n n 


x 5 (t) 


x 5 ( 0 ) cos (nt) + 


x 6 ( 0 ) sin (nt) 
n 


C 5 


( t cot) 


n 


2 


sii^n t) 
2 n 


C 6 t sin (nt) 

2 n 

(4.51) 


xe(t) = x 6 ( 0 ) cos(n£)— nx 5 ( 0 ) sin (nt)-{ 


C 5 t sin (nt) 
2 n 


-C fi 


t cos (nt) sin (nt) 

2 + 2 n 


(4.52) 


Eqs. (4.47)-(4.52) can be solved for the boundary conditions provided to obtain the 
values of constants Ci, C 2 , C 3 , C 4 , C 5 and C 6 , which can then be used to obtain the 
required specific thrust input from Eqs. (4.24)-(4.26). 
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4.3 Optimization of the Control Input Trajectory: 
Co-planar Case 


Considering the maneuvering spaecraft. is at some finite initial co-planar distance 
ro = (xo, yo) relative to the target and there is an initial relative velocity at time 
to = 0 given as v 0 = (v xo , v yo ). 

The thrust is applied along the x— and y —direction as u — (u x , u y ) 

Also, since, we are considering only the co-planar case, the state-space form can be 
expressed as follows: 


Xi 

x 2 

■h 

Xi 


0 

3 n 2 

0 

0 


1 0 
0 0 
0 0 
—2 n 0 


0 

2 n 

1 

0 



Xl 


0 

0 


X2 

+ 

1 

0 


X 3 


0 

0 


X\ 


0 

1 


^X 

u y 


(4.53) 


Here, the interest is to design a fuel optimal control system. 

u = (u x , u y ) is the specific thrust applied by the thrusters onboard the spacecraft 
and due to obvious reasons, it is proportional to its fuel consumption. 

Thus, to minimize the fuel consumption, the performance index to be minimized is 
given as, 


J 


1 

2 


u T Ru dx 


't 0 


(4.54) 


where R G M 3x3 , is a weighing factor such that R T = R. 


The Hamiltonian for the above system with the performance index mentioned in Eq. 
(4.54) is obtained as follows: 

_ l 

"H(T(f), «(£), A(t)) = -u t Ru + \\X 2 + A 2 (3n 2 a;i + 2 nx 4 + u x ) 

+X 3 X 4 + A 4 (- 2 na ;2 + u y ) 


(4.55) 
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u y = -A 4 (4.58) 

Substituting the value of u from Eqs. (4.57) and (4.58) into Eq. (4.55) gives us the 
optimized value of the Hamiltonian as follows: 

77 = \\X2 3 /i X 2 x\ 2n\ 2 X4 3- A 3 X 4 — 2 , 71 X 4 X 2 — “A 2 — 2”^^ (4.59) 

Here, 

1 0 0 

0 0 2 n 

(4.60) 

0 0 1 

2 n 0 0 
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So that, 


e Af = 


0 - — - cos nt 

n n 


4 — 3 cos nt - sin nt 

n 

3 n sin nt cos nt 0 2 sin nt 

—6 nt + 6 sin nt — + - cos nt 1 - sin nt — 3 1 

n n n 

6 n cos nt, — 6 — sin nt 0 4 cos nt — 3 


(4.61) 


As per the necessary conditions for optimality, 


A* = - 


dU 

dx 


= —A A 


T \** 


and, 


x = 


m 

d\ 


n T 


(4.62) 


Now, we have co-state equations independent of state variables expressed as follows: 


Taking, 


1 

* 


0 —3n 2 0 0 


1 - 

* i-H 

a 2 * 


-100 2 n 


A^ 

As* 


0 0 0 0 


As 

1 

* 


0 —2 n -1 0 


K 


(4.63) 


0 —3n 2 0 0 


A = 


-1 

0 


0 

0 


0 2 n 
0 0 


0 —2 n -1 0 


(4.64) 


we obtain 
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e At — 


4 — 3 cos (nt) 
— - sin (nt) 

0 

n(l - cos(nt)) 


—3n sin(nf) 6 nt — 6 sin(nt) 6n cos (nt) — 6n 
cos (nt) ^(cos(nt) — 1) 2sin(nt) 

0 1 0 
—2 sin(nt) 3 1 — - sin (nt) 4 cos (nt) — 3 


(4.65) 


Let the initial values of the co-state variables be, 


A(0) = <<?!, C 2 , C 3 , C 4 ) (4.66) 

Therefore, using, A (t) = e A *A(0), the equations for the co-state variable are obtained 
as follows: 


Ai(t) = Ci(4 — 3 sin nt) + C 2 (—3nsinnt) + C 3 (6nt — 6 sin nt) 

+C 4 (6ncosnf — 6n) (4.67) 


Ci 2C 

A2 (t) =-sin nt + C9 cos nt H-(cos nt — 1) + 2C 4 sin nt 

n n 


(4.68) 


A 3 (t) — C3 


(4.69) 


2C 4 

A 4 (f) = —-(1 — cos nt) — 2 C 2 sin nt, + C 3 (3t -sin nt) + CJAcosnt — 3) (4.70) 

n n 
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Substituting Eqs. (4.68) and (4.70) in Eqs. (4.57) and (4.58), we obtain the optimal 
control trajectories along x- and y- direction as follows: 


c 2C3 

u x (t) = — sin nt — C -2 cos nt - L (cos nt — 1) — 2C 4 sin nt (4.71) 

n n 


2C 4 

u v (t) =--(1 — cos nt) + 2C*2 sin nt — C 3 (3t -sin nt) — C' 4 (4 cos nt — 3) (4.72) 

n n 

With all this information, the solution to the state equation is obtained by, 

e A{t - T) Bu(r) dr (4.73) 

If the initial states, i.e., the initial relative position and velocity of the deputy are 
provided as, 

£(0) = (xi(0) ,x 2 (0), x 3 (0), x 4 (0)) (4.74) 

Multiplying Eqs. (4.61) and (4.74), we obtain: 

4 — 3 cos nt - sin nt 0 - — - cos nt 

n n n 

. 3 n sin nt cos nt 0 2 sin nt 

e At x(0) = 

—6 nt + 6 sin nt — + - cos nt 1 - sin nt — 3 1 

n n n 

6 n cos nt — 6 — sin nt 0 4 cos nt — 3 



x*(t) = e At f(0) + 
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Then, by substituting the value of u from Eqs. (4.71) and (4.72) to the integral part 
of Eq. (4.73) and then integrating it from t = 0 to t and adding Eq. (4.75), the 
solution to the state equation is obtained as follows: 


, , , , 2 . , , w ,4(72 16C 3 

x\(t) = £4(0) — (1 — cos nt) — xi(0)(3cos nt — 4) 4-— 4-— 

n n 2 n 6 

4C 2 13Ci . 11C 4 . 3C 3 2 1 . 4 Ci 

—5- cos nt 4-5- sm nt -— sin nt - 1 4—x 2 (0) sin nt - —t 

n 6 2 n 6 n 2 n n n 2 


6C 4 5Ci 5C 4 5C 2 5 C 3 . 

-t — —1 cos nt 4- 1 cos nt —-—t sm nt - —t sm nt (4.76) 


n 


2 n 2 


n 


2 n 


rr 


x 2 (t) = x 2 (0) cos nt -5- 8C'i — 12C 4 n — 8C\ cos nt — 22C 3 sin nt 

2n 2 \ 

4-12 C 3 nt + 12 C 4 n cos nt — 3 C 2 n sin nt + 10 C 3 nt cos nt — hC 4 nt sin nt 
4-5 C 2 n 2 t cos nt + 10 C 4 n 2 t sin nt\ + 2x 4 (0) sin nt + 3nx 1 (0) sin nt 


(4.77) 


x 3 (t) = x 3 (0) - -x 2 (0)(l - cos nt) - + ^C 4 + - \c 4 t 2 


n 


n° 


rr 


4-6xi (0) (sin nt — 1) — x 4 (0)(3f-sin nt) 4—cos nt 


n 


n° 


28 a 


lic 2 


38C 


A 1 lOo OOO.l 3Cl n UO 2 400,2 

cos nt 4-s— sm nt 4-— sm nt 4- 1 2 - 1 - —t 


6C 2 


28 C 3 


n 2 n 2 n° n n n° 

5C 2 IOC's 5C', 10C 4 , . 

- 1 cos nt -— t cos nt 4-— t sm nt - 1 sm nt (4.78) 


n 


n 


n 


n 
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/X , w , X ^ 6 C 2 28C 3 9C 3 n 

x 4 (t) = x 4 (0)(4cosnt — 3) — 9 C 4 t --—I—— t 2 


n 


rr 


6C 2 


—6naq(0)(l — cos nt) — 2x(0) sin nt H-cos nt 

n 

28C 3 IICi 18C 4 

H-— cos nt -— sin nt H-sin nt 

n 1 n 2 n 

GC 5 c IQ (j 

-10C 4 t cos nt + hC 2 t sin nt H-1 H-1 cos nt H- -t sin nt 

n n n 


(4.79) 


Now, provided we have the initial and final values of the relative position and relative 
velocity of the maneuvering spacecraft, Eqs.(4.76)-(4.79) can be solved to obtain the 
values of the constants C \, C*2, C 3 and C 4 , which can then be substituted into the 
Eqs.(4.71)-(4.72) to obtain the control input trajectories along the respective axes. 


4.4 Coplanar Optimal Relative Control with Sin¬ 
gle Input 


Considering the maneuvering spacecraft is at some finite initial distance r 0 = ( x , y, z) 
relative to the target and there is no relative velocity initially at time to = 0. It is 
required for the maneuvering spacecraft to attain the desired set of distance and 
velocity relative to the target in some specified time tf. 

Since, the only thrust input is along y-axis, we have, u = (0, u y , 0) 

Also, since, we are considering only the co-planar case, the state-space form can be 
interpreted as, 


X\ 


0 

1 

0 

0 


Xi 


0 

X-2 

_ 

3 n 2 

0 

0 

2 n 


x 2 

+ 

0 

x 3 


0 

0 

0 

1 


x 3 


0 

x 4 


0 

—2 n 

0 

0 


x 4 


1 


(4.80) 
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subjected to performance index, 


J = 



so that, the Hamiltonian is obtained as follows: 


(4.81) 


H = -Uy + Xix 2 + A 2 (3n 2 a;i + 2 nx 4 ) + A 3 X 4 + A 4 (—2 nx 2 + u y ) (4.82) 


Here, 


A 


0 10 0 
3n 2 0 0 2 n 

0 0 0 1 
0 — 2 n 0 0 


So that, 


„At 


4 — 3 cos nt - sin nt 

n 

3 n sin nt cos nt 

- 6 nt + 6 sin nt — + 

n n 

6 n cos nt — 6 — sin nt 


0 

0 

1 

0 


22 j . 

-cos nt 

n n 

2 sin nt 

4 sin nt _ 

n 

4 cos nt — 3 


(4.83) 


(4.84) 


As per the necessary conditions for optimality, 


A* = 


dU 

dx 


= -A 1 A 


T \ * 


Now, we have co-state equations independent of state variables obtained as follows: 


1 

* 


1 

0 

1 

GO 

3 

to 

0 

0 


1 - 

* 1—1 

x 2 * 


-100 2 n 



a 3 * 


0 0 0 0 


A 3 

1 

* 


0 -2 n -1 0 


_A1_ 


(4.85) 


Taking, 
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0 —3 n 2 0 0 


A = 


-1 0 

0 0 


0 2 n 
0 0 


0 -2 n -1 0 


Assuming initial values for X T (0) = [Ai(0), A2(0), As(0), A4(0)] = C = [C\, 
and as, 


A *(t) = e At C 


The co-state equations are obtained as follows: 


Ai(t) = C\ (3 cos rd — 4) — C' 3 (6 sin nt, — 6 nt) — C^Gn cos nt — 6) 

—3C 2 n sin nt 


\ 2 (t) 



2 cos nt 


n 


— C 2 cos nt + C 4 sin nt — C\ — sin nt 

n 


(4.86) 


■ 2 , CM], 


(4.87) 


(4.88) 


(4.89) 


A 3 (t) = -C :i 


(4.90) 
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4 sin Tit 

A 4 (t) = C 3 [ 3t -I — C 4 (4 cos nt — 3) — 2C 2 sinnf 


n 


.2 2 cos nt 

-Li - 


n n 


(4.91) 


The stationarity (necessary condition for optimality) is as follows: 


dn _ _ 

— 0 — n t + X l 


which results in, 


„ * \ * 

Uy — — A 4 


Substituting the value of A 4 from Eq. (4.91) into above equation, we obtain the 
expression for optimal control along y-axis as follows: 


4 sin nt \ „ ___ „ / 2 2 cos nt 


n 


u*,. = — C 3 ( 3 1 -] + C 4 (4 cos nt — 3) + 2C 2 sin nt + Ci I — 


) 


\ 


n n 


(4.92) 


The states are obtained as follows: 


x(t) = e At x{t}) + J e A(t ~ T) Bu(r)dT 


(4.93) 


If we consider the initial values of the states of the state-space model as, 


£(0) = [:£i(0), x 2 (0), :r 3 (0), £ 4 (0 )] t 


(4.94) 
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which gives, 


e At x{0) 


4 — 3 cos nt 

- sin nt 

n 

0 

2 2 + 
-cos nt 

n n 


®i( 0 ) 

3 n sin nt 

cos nt 

0 

2 sin nt 


2:2(0) 

—6 nt + 6 sin nt 

—2 | 2 cos nt 

n n 

1 

4 sin nt 

n 


2:3(0) 

6 n cos nt — 6 

— sin nt 

0 

4 cos nt — 3 


2:4(0) 


(4.95) 


Hence, by solving Eq. (4.93) along with the initial values of position and velocity give 
the solution to the optimal states and the along the track optimal input trajectory 
can be obtained. 


4.5 Optimal Low-Thrust Constrained Optimiza¬ 
tion 


Considering the maneuvering spaecraft- is at some finite initial co-planar distance 
f*o = (xo, Do) relative to the target and there is an initial relative velocity at time 
t 0 = 0 given as v 0 = (v xo , v yo ). 

The thrust is applied along the x— and y —direction as u— (u x , u y ) 

Also, since, we are considering only the co-planar case, the state-space form can be 
expressed as follows: 


Xi 

X2 

■h 

CC4 


0 

3 n 2 

0 

0 


1 0 0 
0 0 2 n 

0 0 1 
—2 n 0 0 


2:1 


0 

0 

2:2 

+ 

1 

0 

2:3 
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0 

£4 


0 

1 


u x 

Uy 


(4.96) 


Here, the interest is to design a fuel optimal control system. 

u = (u x , u y ) is the specihc thrust applied by the thrusters onboard the spacecraft, 
and is constrained within the limits: 
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u(t) | < r t n 2 


(4.97) 


Thus, to minimize the fuel consumption, the performance index to be minimized is 
given as, 


J 


1 

2 


u T Ru dx 


'to 


(4.98) 


where R € M 3x3 , is a weighing factor such that R T 


R. The optimal control is: 


u*(t) = —sat{g*(f)} (4.99) 

where, 

q* = R~ l B T \(t) (4.100) 


Now, proceeding in the same way as in the case of unconstrained problem in Section 
(4.3), we have, the co-state equation obtained as: 

Ai(t) = C*i (3 cos nt — 4) — Cs(6 sin nt — 6 nt) — C±{Qn cos nt — 6) 

—3C 2 nsinni (4.101) 


A 2 (t) 



2 cos nt 


n 


C 2 cos nt + C 4 sin nt — C\ — sin nt 

n 


(4.102) 


A 3 (t) = -C 3 


(4.103) 
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4 sin Tit 

A 4 (t) = C 3 I 3t - ) — (^4(4 cos nt — 3) — 2U 2 sinnt 


n 


. 2 2 cos nt 

-Ci - 


n n 


(4.104) 


Using Eq. (4.100), we obtain: 


2 1 

qUt) = C 2 cos (nt) -U 3 cos (nt) + 2 C' 4 sin(nt)- C\ sin (nt) (4.105) 

n n 


2 4 

q 2 (t) = — CUl — cos (nt)) + 6*4 (4 cos (ni) — 3) — 2U 2 sin(nt) + C 3 (3t -sin(nt)) 

n n 

(4.106) 


Hence, from Eq. (4.99), we obtain: 


u*(t) = { 


r t n 2 , 


—r t n , 


if q*(t) < —r t n 2 . 
if q*(t) > r t n 2 . 
if \q*(t)\ < r t n 2 . 


(4.107) 


The above system of equations can be numerically solved to obtain the optimal 
control input trajectories. 



Chapter 5 


Numerical Analyses 


The following chapter presents the numerical analyses of the optimization process 
discussed in Chapter 4. Various cases have been considered here for the optimization 
such as the coplanar rendezvous of a maneuvering spacecraft with the ISS and 
rendezvous for a more general case considering the out-of-plane motion also. 
Furthermore, the comparison of the optimal control input trajectory as obtained 
with the HCW model has made with that obtained using the non-linear model. Since, 
the inclusion of the oblateness effect adds a finite error in the control input obtained 
using HCW model, a comparison has been made between the control input trajectory 
obtained with and without considering the obtlateness effects of the Earth. 

For numerical analyses, MATLAB IM inbuilt solver bvp4c has been used. 
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5.1 Optimal Control Input Trajectory for the Copla- 
nar Case 

The following section shows the results of the numerical simulation carried out for 
the coplanar motion, he., considering motion only along radial (x) and along track 
(y) direction only. 

5.1.1 Only In-Track Control Input Application 

Case 1 

Rendezvous Problem: Consider a simple case, when a spacecraft at a radial distance 
of 100 Km from International Space Station (ISS) has to rendezvous with it within a 
period of {it/ n) seconds, where ISS is located at an approximated altitude of 400 
Km above Earth surface, having orbital frequency n = 1.1313658 x 10 -3 . 

Considering only co-planar case, we have following data, 

x{t = 0) = xi(0) = 100 Km 

v x {t = 0) = £ 2 ( 0 ) = 0 Km/s 

y{t = 0) = £ 3 ( 0 ) = 0 Km 

v y (t = 0) = 2 : 4 ( 0 ) = 0 Km/s, and 

tf = 7T /n s 

Figure (5.1) shows the plot of required control input with the time of the rendezvous 
operation, which is of the order of 10 -3 Km/s 2 . Figure (5.2) shows the respective 
trajectories of the relative distances along x and y direction, and it can be observed 
that at the specified time of t — tt/u seconds, both the relative distance components 
are settcling to zero values, showing the successful rendezvous. Same result can be 
predicted from the figure (5.3) which is showing the respective relative velocities 
trajectories with time. Finally, figure (5.4) shows the coplanar geometry of the 









Time (s) 
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Case 2 

Rendezvous Problem: Continuing with the case 1, but now considering that the 

maneuvering spacecraft is having some initial velocities also, 

x(t = 0) = xi(0) = 100 Km 

v x (t = 0) = £ 2 ( 0 ) = 0.02 Km/s 

y{t = 0) = £ 3 ( 0 ) = 0 Km 

v y (t = 0) = £ 4 ( 0 ) = 0.02 Km/s, and 

tf = n/n s 

When the maneuvering spacecraft is having a finite initial relative velocity value 
relative to the target, the required control input increases as compared to the Case 1, 
which is also obvious, as shown in figure (5.5). Figure (5.6) and figure (5.7) shows the 
trajectories of the relative distances and relative velocities respectively. Figure (5.8) 
shows the geometry of the orbit of the maneuvering spacecraft around the target. 



Specific Thrust (Km/s 2 ) xIO 


Figure 5.5: Trajectory of Control Input Along y-Direction: Case 2 
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5.1.2 Control Input Applied Along Radial and Along-Track 
Direction 

Case 3 

Rendezvous Problem: Considering a case when the maneuvering spacecraft has to 
rendezvous with the ISS, located at a radial distance of 50 Km and in-track distance 
of 50 Km. In this case, it is assumed that the maneuvering spacecraft is initially 
stationary relative to ISS. 

Hence, considering only co-planar case, we have following data, 

x(t = 0) = 2 q( 0 ) = 50 Km 

v x (t = 0) = 2 : 2 ( 0 ) = 0 Km/s 

y(t = 0) = 2 : 3 ( 0 ) = 50 Km 

v y (t = 0) = 2 : 4 ( 0 ) = 0 Km/s, and 

tf = 7T /n s 

Figures (5.9) and (5.10) shows the trajectory of the control input along x and y- 
direct.ion, respectively. It can be seen that the order of the required control input 
is 10 -4 Km/s 2 . Figure (5.11) shows the trajectories of the relative distances along 
x and y-axis, respectively, while figure (5.12) shows the relative distance of the 
maneuvering spacecraft from the target. Figure (5.13) shows the trajectories of the 
relative velocities along x and y-axis. Figure (5.14) shows the relative orbit of the 
maneuvering spacecraft around the target spacecraft. 
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Case 4 

Rendezvous Problem: Considering a case when the maneuvering spacecraft has to 
rendezvous with the ISS, located at a radial distance of 50 Km and in-track distance 
of 50 Km. In this case, it is assumed that the maneuvering spacecraft is having an 
initial relative velocity of 0.025 Km/s, along both radial and tangential directions 
relative to ISS. 

Hence, for this case, we have following data, 

x(t = 0) = a;i(0) = 50 Km 

v x (t = 0) = £ 2 ( 0 ) = 0.025 Km/s 

y(t = 0) = £ 3 ( 0 ) = 50 Km 

v y (t = 0) = £ 4 ( 0 ) = 0.025 Km/s, and 

tf = tt/ n s 

Figures (5.15) and (5.16) shows the trajectory of the control input along x and 
y-direction, respectively. It can be seen that the order of the required control input 
is 10 -4 Km/s 2 . Figure (5.17) shows the trajectories of the relative distances along 
x and y-axis, respectively, while figure (5.18) shows the relative distance of the 
maneuvering spacecraft from the target. Figure (5.19) shows the trajectories of the 
relative velocities along x and y-axis. Figure (5.20) shows the relative orbit of the 
maneuvering spacecraft around the target spacecraft. 
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Case 5 

Stationkeeping Problem: Considering a case when the maneuvering spacecraft has to 
maintain a distance of 10 Km relative to the ISS along both radial and tangential 
direction, which is initially located at a radial distance of 50 Km and in-track distance 
of 50 Km. In this case, it is assumed that the maneuvering spacecraft is having an 
initial relative velocity of 0.025 Km/s, along both radial and tangential directions 
relative to ISS. 

Hence, for this case, we have following data, 

x(t = 0) = xi(0) = 50 Km 

v x (t = 0) = £ 2 ( 0 ) = 0.025 Km/s 

y(t = 0) = £ 3 ( 0 ) = 50 Km 

v y (t = 0) = £ 4 ( 0 ) = 0.025 Km/s, and 

tf — n/n s 

Figures (5.21) and (5.22) shows the trajectory of the control input along x and 
y-direction, respectively. It can be seen that the order of the required control input is 
10 -4 Km/s 2 . Figure (5.23) shows the trajectories of the relative distances along x and 
y-axis, respectively, while figure (5.24) shows the relative distance of the maneuvering 
spacecraft from the target. It can be seen that at the specified time, the maneuvering 
spacecraft is at the specified distance from the target, however, meanwhile it does 
come closer than the specified relative distance. Figure (5.25) shows the trajectories 
of the relative velocities along x and y-axis. Figure (5.26) shows the relative orbit of 
the maneuvering spacecraft around the target spacecraft. 
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5.2 Comparison of the Optimal Control Input Tra¬ 
jectories as Obtained using HCW and Non- 
Linear Equations for the Coplanar Case 


To examine the limit of compatibility of the HCW equations, a comparison with 
the optimal input trajectory obtained using the non-linear equations of motion is 
necessary. The following cases compares the results as obtained using the HCW 
equations and the non-linear equations. 

5.2.1 Comparison for the Case of Rendezvous 

Considering the case of rendezvous of the maneuvering spacecraft with the ISS, 
located at an altitude of 400 Km above Earth’s surface, the relative distance of 
the maneuvering spacecraft is 50 Km along both radial and tangential direction. 
Comparison for the results has been made for different mission timings, t = n/A:n 
and t = 7t/2 n seconds. 

Case 6 

If the maneuvering spacecraft is having an initial relative velocity of 0.025 Km/s along 
both radial and tangential directions, and it is required to complete the rendezvous 
in 7t/4 n seconds, the comparison of applied control inputs along both x- and y- 
direction are given here. 

Figure (5.27) shows the trajectory of the relative distances. It can be seen that the 
relative distance becomes 70.71 Km at the specified distance. It can be seen that 
the solutions for the trajectories are coinciding and perfectly in agreement. Figures 
(5.28) and (5.29) shows the comparison of the control input trajectories. In this case, 
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some deviation can be observed from the non-linear solution. Figures (5.30) and 
(5.31) shows the respective deviation of the control inputs along x and y-axis, which 
is of the order of 10 -5 Km/s 2 . Figures (5.32) and (5.33) shows the comparison of 
the relative velocity trajectories along x and y-axis. 



Figure 5.27: Trajectory of Relative Distances: Case 6 





Figure 5.28: Trajectory of Control Inputs Along x-Direction: Case 6 



Figure 5.29: Trajectory of Control Inputs Along y-Direction: Case 6 
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- 0.12 - 0.1 - 0.08 - 0.06 - 0.04 - 0.02 0 0.02 

v x (Km/s) 

Figure 5.32: Comparison of Relative Velocity Component Along x-Direction: Case 6 



- 0.08 - 0.06 - 0.04 - 0.02 

v y (Km/s) 


Figure 5.33: Comparison of Relative Velocity Component Along y-Direction: Case 6 
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Case 7 

If the maneuvering spacecraft is having an initial relative velocity of 0.025 Km/s 
along both radial and tangential directions, but now the time of rendezvous has been 
increased to n/2n seconds, the comparison of applied control inputs along both x- 
and y- direction are given here. 

Figure (5.34) shows the trajectory of the relative distances. It can be seen that the 
relative distance becomes 70.71 Km at the specified distance. It can be seen that, 
unlike the previous case, in this case the solutions for the trajectories are showing 
slight variation. Figures (5.35) and (5.36) shows the comparison of the control input 
trajectories. In this case, some deviation can be observed from the non-linear solution. 
Figures (5.37) and (5.38) shows the respective deviation of the control inputs along 
x and y-axis, which is of the order of 10~ 5 Km/s 2 . Figures (5.39) and (5.40) shows 
the comparison of the relative velocity trajectories along x and y-axis. 



Figure 5.34: Trajectory of Relative Distances: Case 7 
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Figure 5.35: Trajectory of Control Inputs Along x-Direction: Case 7 
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Figure 5.36: Trajectory of Control Inputs Along y-Direction: Case 7 
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5.3 Comparison of the Results for Coplanar Low- 
Thrust Rendezvous Optimization with HCW 
and Non-Linear Model 

This section tabulates the the respective maximum errors in relative velocity, relative 
distance and required control input found by comparing the simulation results of 
HCW and non-linear model, when compared with the different sets of specified initial 
conditions and time of operation for the coplanar rendezvous obtained by studying 
different cases as in previous section. 

Table 5.1 shows the results of the different cases of coplanar rendezvous of maneuvering 
spacecraft located initially at a relative distance of 70.71 Km from the ISS, for a 
specified operation time of {it/ n) = 2776.8 seconds and provided with different initial 
velocities. 


Wo = Vyo (Krn/s) 

^^max (Km/s) 

A r max (Km) 

A Umax (Km/s 2 ) 

0 

0.0025 

4.16 

0.0001463 

0.025 

0.004136 

7.4 

0.0002092 

0.05 

0.00644 

10.19 

0.0002714 

0.1 

0.01146 

16.2 

0.0003951 

0.2 

0.02728 

27.8 

0.000638 


Table 5.1: Coplanar Low-Thrust Rendezvous Optimization with HCW and Non-Linear 
Model: Varying the Initial Relative Velocitites 


It can be deduced from the data given in table 5.1 that as the initial relative velocity 
of the maneuvering spacecraft increases, there is increase in the error in the relative 
velocity, relative distance and the optimal control input. 
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Table 5.2 shows the results of the different cases of coplanar rendezvous of maneuvering 
spacecraft with ISS in a specified operation time of (n/n) = 2776.8 seconds, with 
initial relative velocity of 0.025 Krn/s and with different sets of relative distances 
from ISS. 


x 0 = Vo (Km) 

AtWr (Krn/s) 

A r max (Km) 

A u m a X (Km/s 2 ) 

50 

0.0002092 

7.15 

0.004136 

100 

0.01282 

11.26 

0.0003481 

200 

0.04342 

19.4 

0.000609 

300 

0.0915 

27.5 

0.000845 

500 

0.2187 

42.7 

0.001251 


Table 5.2: Coplanar Low-Thrust Rendezvous Optimization with HCW and Non-Linear 
Model: Varying the Initial Relative Distances 


From the table 5.2, it can be seen that as the initial relative distance of the ma¬ 
neuvering spacecraft relative to the target spacecraft increases, there is a drastic 
increment in the errors in the relative velocity and it becomes impractical to use the 
HCW model for optimization beyond 100 Km of relative separation, since, the error 
becomes more than 50% of the initial relative velocity. 
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Table 5.3 shows the results of the different cases of coplanar rendezvous of maneuvering 
spacecraft with ISS. The maneuvering spacecraft is initially located at a relative 
distance of 70.71 Km from the ISS and having an initial relative velocity of 0.025 
Krn/s relative to ISS. The simulations were performed for different sets of specified 
time of operation. 


t (s) 

Av max (Km/s) 

A r max (Km) 

A u m a X (Km/s 2 ) 

7r/10n 

0.0003819 

0.020135 

0.0024666 

7t/4 n 

0.0004953 

0.1934 

0.00049 

n/3n 

0.000580 

0.576102 

0.000353 

7t/2 n 

0.000868 

2.40901 

0.000281 


Table 5.3: Coplanar Low-Thrust Rendezvous Optimization with HCW and Non-Linear 
Model: Varying the Specified Time of Operation 


From the table 5.3, it can be seen that as the specified time of operation increases, 
error in the relative velocities and relative distance increases, however, error in control 
input as predicted by the HCW model decreases. This can be contributed to the 
fact that longer time of operation lowers the continuous thrust requirements and so 
as the errors in the HCW model. 
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5.4 Comparison of the HCW and Non-Linear Mod¬ 
els for the Case of Stationkeeping 


Case 8 

Considering the case when the maneuvering spacecraft has to maintain a distance of 
10 Km along both tangential and radial direction i.e., a relative distance of 14.14 Km 
with the target (here, ISS), which is in a circular orbit of radius 6778.14 Km. The 
plots here shows the results of the simulation carried out for operation time from 
t — 50 s to t — ( ir/n ) s and minimum relative distance reached has been plotted in 
each case. The maneuvering spacecraft is initially at relative distance of 70.71 Km 
and having relative velocities 0.025 Km/s along both radial and tangential directions 
relative to the ISS. 

It can be seen from the figure (5.41) that as per the optimization using the non-linear 
equations of relative motion, at and beyond time of operation 1600 s, the maneuvering 
spacecraft crosses the specified station-keeping distance, thereby, making the mission 
vulnerable to collision. However, the HCW model predicts the same at 1760 s. 
Furthermore, if we examine the plots for maximum specific thrust requirements for 
the case, it can be observed that the specific thrust requirement decreases with the 
increase in the specified mission timing. 
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Case 9 

Continuing with the above case 8, with slight change in the initial conditions for 
the maneuvering spacecraft, increasing the initial relative velocities (along both 
radial and tangential direction) to 0.05 Km/s, the simulation has been performed for 
different operation timings and results have been presented in the plots. 

It can be observed from the figure (5.43) that in this case, as predicted by the 
non-linear relative equations model that the maneuvering spacecraft crosses the 
specified station-keeping distance for operation timings beyond 1550 s. However, 
the results varies when the results of the HCW simulation are taken into account, 
according to which, the same limit is 1690 s. 

It can be further observed that the limit on time of operation decreases as the initial 
relative velocities increases. To further study this, simulation has been repeated for 
further increased initial relative velocity in Case 10. 



min 



Figure 5.43: Comparison of the Minimum Relative Distances Reached: Case 9 



t (S) 


Figure 5.44: Comparison of the Peak Thrust Requirements: Case 9 
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Case 10 

This case is similar to the one presented in the Case 8, however, the initial relative 
velocities have been increased to 0.1 Km/s. The simulation has been performed for 
different operation timings and results have been presented in the plots. 

It can be observed from the figure (5.45) that in this case, as predicted by the 
non-linear relative equations model that the maneuvering spacecraft crosses the 
specified station-keeping distance for operation timings beyond 1490 s. However, 
the results varies when the results of the HCW simulation are taken into account, 
according to which, the same limit is 1620 s. 




Figure 5.46: Comparison of the Peak Thrust Requirements: Case 10 
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5.4.1 Summary of the Results: Optimal Time for Station¬ 
keeping 


From the studied cases (8-10) for different initial velocities specified, it can be observed 
that the required control effort decreases with time. So, one way of obtaining the 
optimal time for operation is to opt for a longer duration. However, it is evident 
from the plots that beyond a certain time of operation, the relative distance between 
maneuvering spacecraft and the target becomes less than the specified nominal 
distance of 14.14 Km. Hence, the optimal time of operation is restricted by the 
time-limit when the actual relative distance between the two spacecraft becomes less 
than the specified nominal distance. Table 5.4 enlists the optimal time of operation 
for different initial velocities obtained by the above method. 


u(0) 

toptimal (s) 

(Km/s) 

HCW 

Non-Linear 

0.025 

1760 

1600 

0.05 

1690 

1550 

0.1 

1620 

1490 


Table 5.4: Stationkeeping: Optimal Time for Different Initial Speeds 
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5.5 Optimal Control Input Trajectory for the Gen¬ 
eral Case 


Whereas the previous sections deals with the coplanar relative motion in the orbit, 
this section onwards will deal with the cases that taken in account the out-of-plane 
motion i.e ., the motion along 2 -axis also. 

Case 11 

Rendezvous Problem: Consider that the maneuvering spacecraft has to rendezvous 

with the ISS, which is in a circular orbit at an altitude of 400 Km above Earth’s 

surface. For this particular case, we have the following data, 

x(t = 0) = xi(0) = 50 Km 

v x (t = 0) = £ 2 ( 0 ) = 0 Km/s 

y(t = 0) = £ 3 ( 0 ) = 50 Km 

v y (t = 0) = £ 4 ( 0 ) = 0 Km/s 

z(t = 0) = £ 5 ( 0 ) = 50 Km 

tq(0) = a; 6 (0) = 0 Km/s 

and tf = 7 r /n s 

Figure (5.47) shows the trajectory of the control input which is of the order of 10 -4 
Km/s 2 . Figure (5.48) shows the respective relative distance trajectories along x, 
y and 2 axis, each starting from an initial separation of 50 Km and reaching the 
required zero distance in the specified time. Figure (5.49) shows the trajectory of 
the relative distances along the respective x,y and 2 axis. The initial and final values 
of the relative distances can be noticed, which are as per the specified initial relative 
distances. Figure (5.50) shows the geometry of the maneuvering spacecraft orbit. 
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Case 12 

Rendezvous Problem: Consider that the maneuvering spacecraft has to rendezvous 

with the ISS, which is in a circular orbit at an altitude of 400 Km above Earth’s 

surface. For this particular case, it has finite initial relative velocity relative to the 

ISS, so that, we have the following data, 

x(t = 0) = 2 q( 0 ) = 50 Km 

v x (t = 0) = 2 : 2 ( 0 ) = 0.025 Krn/s 

y(t = 0) = 2 : 3 ( 0 ) = 50 Km 

v y (t = 0) = 2 : 4 ( 0 ) = 0.025 Krn/s 

z(t = 0) = 2 : 5 ( 0 ) = 50 Km 

u 2 (0) = 2 : 6 ( 0 ) = 0.025 Km/s 

and tf = 7 T /n s 

Figure (5.51) shows the trajectory of the control input which is of the order of 10 -4 
Km/s 2 . Figure (5.52) shows the respective relative distance trajectories along x, 
y and z axis, each starting from an initial separation of 50 Km and reaching the 
required zero distance in the specified time. Figure (5.53) shows the trajectory of 
the relative distances along the respective 2 :,y and z axis. The initial and final values 
of the relative distances can be noticed, which are as per the specified initial relative 
distances. Figure (5.54) shows the geometry of the maneuvering spacecraft orbit. 



Relative Distance (Km) i Control Input (Km/s 


-1 

1.5 

-2 

2.5 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Fraction of Specified Time 

Figure 5.51: Trajectory of Control Input: Case 12 
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Figure 5.52: Trajectory of Relative Distance: Case 12 






z (« m ) Relative Velocity 








0.7 


0.8 




Figure 5.54: Geometry of the Maneuvering Spacecraft Orbit: Case 12 
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5.6 Comparison of Numerical Solution and the 
Closed-Form Solution 

To check the accuracy of the closed-form solution obtained in chapter 4 for obtaining 
the specific thrust requirement for the fuel optimal case, it has been compared with 
the numerical solution obtained using MATLAB bvp4c scheme numerically. 

Case 13 

Rendezvous Problem: For the case of rendezvous with the target which is in a circular 

Earth orbit of radius 6900 Km. For this particular case, we have the following data, 

x(t = 0) = aq(0) = 50 Km 

v x (t = 0) = £ 2 ( 0 ) = 0 Km/s 

y(t = 0) = £ 3 ( 0 ) = 50 Km 

v y (t = 0) = £ 4 ( 0 ) = 0 Km/s 

z(t = 0) = £ 5 ( 0 ) = 50 Km 

u z (0) = a; 6 (0) = 0 Km/s 

and tf = tt/u s 


Substituting above data in state-solutions obtained in chapter 5, viz. equations (4.80) 
to (4.85), and solving for the values of constants C\, C 2 , C 3 , C 4 , C 5 and C 6 , we get, 
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Ci = 4.77 x 10 -7 (5.1) 

C 2 = 1.42 x 10~ 4 (5.2) 

C 3 = -4.79 x 10" 8 (5.3) 

C 4 = 1.83 x 10~ 4 (5.4) 

C 5 = 4.25 x 10" 8 (5.5) 

C 6 = 1.09 x 10“ 36 (5.6) 

Substituting these values in equations (4.57), (4.58) and (4.59), we get, 

u x (t) = 0.00006544 sin(0.00110152t) 

-0.00005595 cos(0.001101521) - 0.00008700 (5.7) 

u y (t) = 0.000000141 + 0.00013089 cos(0.00110152t) 

+0.00011190 sin(0.001101521) - 0.00031485 (5.8) 

u z (t) = 0.00003862 sin(0.00110152 1) 

- (1.089 x 10” 36 ) cos(0.00110152 1) (5.9) 
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Also, the solution to the states can be obtained from equations (4.80) to (4.85), 


X!(t) = 653.82 sin(O.OOllt) -47.64 cos(0.0011 1) - 0.5716 1 
-0.1485 1 cos(O.OOllt) - 0.1269 1 sin(0.0011 1) + 0.00013 1 2 + 97.64 (5.10) 

x 2 (t) = 0.00026 t + 0.5716 cos(0.0011 t) - 0.0745 sin(O.OOllt) 

-0.00014 1 cos(O.OOllt) + 0.00016 1 sin(0.0011 1) - 0.5716 (5.11) 

x 3 (t) = 1469.46 cos(0.0011 1) — 0.0033 1 + 233.64 sin(O.OOllt) 

—0.25391 cos(0.0011 1) + 0.2970 1 sin(0.0011 1) + 0.00047 1 2 

-0.000000072 1 3 - 1419.46 (5.12) 


x 4 (t) = 0.00094 1 + 0.00337 cos(O.OOllt) - 1.3215 sin(O.OOllt) 
+0.00032 1 cos(O.OOllt) + 0.00027 1 sin(O.OOllt) - 0.00000021 1 2 - 0.00337 (5.13) 

x 5 (t) = 50.0 cos(O.OOllt) + 15.9154 sin(O.OOllt) - 0.0175 1 cos(O.OOllt) 

- (4.9444 • 10 -34 ) t sin(0.0011 1) (5.14) 


x 6 (t) = 0.000019 1 sin(O.OOllt) - (5.4464 • 10“ 37 ) t cos(O.OOllt) 

-0.0550 sin(O.OOllt) 


(5.15) 



(Km/s 
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The plots below shows the comparison of the results obtained using the closed-form 
solution and the numerical solution. Figures (5.55) to (5.57) shows the comparison 
of the control inputs along respective x, y and z-axis. It can be seen that both the 
solutions are in close agreement. Figures (5.58) to (5.60) shows the comparison of 
the relative distances along x, y and z-axis. It can be seen that both the solutions 
of the relative distances are in close agreement. Figures (5.61) to (5.63) shows the 
comparison of the relative velocities along x, y and z-axis. It can be seen that both 
the solutions are in close agreement. 



Figure 5.55: Comparison of the Trajectory of Control Input Along x-direction: Case 13 





Figure 5.57: Comparison of the Trajectory of Control Input Along z-direction: Case 13 







Figure 5.59: Comparison of the Relative Distance Trajectory Along y-direction: Case 13 
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Figure 5.60: Comparison of the Relative Distance Trajectory Along z-direction: Case 13 



Figure 5.61: Comparison of the Relative Velocity Trajectory Along x-direction: Case 13 
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Case 14 

The following instance has been taken from Marinescu(1976). For the case of 

rendezvous with the target which is in a circular Earth orbit of radius 6678 Km. For 

this particular case, we have the following data, 

x(t = 0) = £i(0) = 2 Km 

v x {t = 0) = £ 2 ( 0 ) = —0.008 Km/s 

y{t = 0) = £ 3 ( 0 ) = —9 Km 

v y (t = 0) = £ 4 ( 0 ) = 0.04 Km/s 

z(t = 0) = £ 5 ( 0 ) = 0.9 Km 

rq(0) = iCg(0) = —0.004 Km/s 

and tf = 2700s 

Substituting values from equations (5.1) and (5.2) in state-solutions obtained in 
chapter 5, viz. equations (4.60) to (4.65), and solving for the values of constants C \, 
C 2 , C 3 , C 4 , C 5 and C 6 , we get, 


Cl = 1.80 x 10' 7 

(5.16) 

C 2 = 4.02 x 10~ 5 

(5.17) 

C 3 = -1.13 x 10' 8 

(5.18) 

C 4 = 7.44 x 10' 5 

(5.19) 

C 5 = 8.86 x 10“ 10 

(5.20) 

C 6 = -2.97 x 10" 6 

(5.21) 
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Substituting these values in equations (4.57), (4.58) and (4.59), we get, 

u x (t) = 7.309 x 10“ 6 sin(1.156 x 10~ 3 1) 

-2.055 x 10~ 5 cos(1.156 x 10“ 3 t) - 1.970 x 10“ 5 (5.22) 

Uy(t) = 3.418 x 10~ 8 t + 1.462 x 10~ 5 cos(l,156 x 10" 3 t) 

+4.110 x 10" 5 sin(1.156 x 10" 3 t) - 8.908 x 10“ 5 (5.23) 

u z (t) = 2.979 x 10~ 6 cos(l.l56 x 10” 3 t) 

+7.665 x 10" 7 sin(l,156 x 10” 3 t) (5.24) 

Also, the solution to the states can be obtained from equations (4.80) to (4.85), 

aq (t) = 139.85 sin(O.OOllt) - 77.68 cos(O.OOllt) - 0.15 1 
-0.015 1 cos(0.0011 1) - 0.044 1 sin(0.0011 1) + 0.000029 1 2 + 79.68 (5.25) 

x 2 (t) = 0.000059 1 + 0.14 cos(0.0011 1) + 0.045 sin(0.0011 1) 

-0.000051 1 cos(0.0011 1) + 0.000018 1 sin(0.00m) - 0.15 (5.26) 

x 3 (t) = 296.08 cos(O.OOllt) - 0.104 1 + 201.43 sin(O.OOllt) 

—0.08 1 cos(0.0011 1) + 0.031 1 sin(0.0011 1) 

+0.00013 1 2 - 0.000000017 1 3 - 305.08 (5.27) 
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x 4 (t) = 0.000261 + 0.14 cos(0.00111) - 0.31 sin(0.00111) 

+0.0000361 cos(O.OOlll) + 0.00011 sin(O.OOlll) 

-0.0000000511 2 - 0.104 (5.28) 

x 5 (t) = 0.9 cos(O.OOlll) - 3.17 sin(O.OOlll) 

-0.000331 cos(0.00111) + 0.00121 sin(0.00111) (5.29) 

x 6 (t) = 0.00024 sin(O.OOlll) - 0.004 cos(O.OOlll) 

+0.00000141 cos(O.OOlll) + 0.000000381 sin(O.OOlll) (5.30) 


The plots below shows the comparison of the results obtained using the closed-form 
solution and the numerical solution. Figures (5.64) to (5.66) shows the comparison 
of the control inputs along respective x, y and z-axis. It can be seen that both the 
solutions are in close agreement. Figures (5.67) to (5.69) shows the comparison of 
the relative distances along x, y and z-axis. It can be seen that both the solutions 
of the relative distances are in close agreement. Figures (5.70) to (5.72) shows the 
comparison of the relative velocities along x, y and z-axis. It can be seen that both 
the solutions are in close agreement. 




Figure 5.64: Comparison of the Trajectory of Control Input Along x-direction: Case 14 
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Figure 5.65: Comparison of the Trajectory of Control Input Along y-direction: Case 14 
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Figure 5.66: Comparison of the Trajectory of Control Input Along z-direction: Case 14 



Figure 5.67: Comparison of the Relative Distance Trajectory Along x-direction: Case 14 







Figure 5.68: Comparison of the Relative Distance Trajectory Along y-direction: Case 14 



Figure 5.69: Comparison of the Relative Distance Trajectory Along z-direction: Case 14 
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Figure 5.70: Comparison of the Relative Velocity Trajectory Along x-direction: Case 14 



Figure 5.71: Comparison of the Relative Velocity Trajectory Along y-direction: Case 14 
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5.7 Comparison of Results Obtained with Non- 
Linear Equations and HCW Equations 


To compare the applicability of the HCW model with the Non-Linear equations of 
relative motion in the orbit, a number of cases are considered for different specified 
time of operation. In the studied cases, the maneuvering spacecraft has to ren¬ 
dezvous with the ISS (a = (6378.14 + 400) Km). It is assumed that the components 
of initial and final relative velocity of the maneuvering spacecraft are zero relative to 
target (ISS) as well as the final relative position. 

For these particular case, we have the following data, 

x(t = 0) = aq(0) = 100 Km 

v x {t = 0) = 2 : 2 ( 0 ) = 0 Km/s 

y(t = 0) = 2 : 3 ( 0 ) = 100 Km 

v y (t = 0) = 2 : 4 ( 0 ) = 0 Km/s 

z(t = 0) = 2 : 5 ( 0 ) = 100 Km 

rq(0) = 2 +( 0 ) = 0 Km/s 
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Case 15: t = n/An 

For specified time of t — n/An = 694.2s, figure (5.73) shows the error in control input 
trajectory. It can be seen that the maximum error in control input is 1.6 x 10” 3 
m/s 2 along the x-axis and 6.68 x 10 -3 m/s 2 along both y and z-axis. Figure (5.74) 
shows the error in relative position, the maximum errors being 18.23 m along radial, 
7.33 m along the track and 11.96 m along z-axis. Figure (5.75) shows the error in 
relative velocity. The maximum error in relative velocity in this case is 0.092 m/s 
along radial, 0.07 m/s along track and 0.097 m/s along the z-axis. 



Figure 5.73: Error in the Control Input: Case 15 




Error in Relative Velocity (m/s) 



Figure 5.74: Error in the Relative Position: Case 15 
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Case 16: t = 7r/3 n 

For specified time of t — ir/3n = 925.6s, figure (5.76) shows the error in control input 
trajectory. It can be seen that the maximum error in control input is 2.4 x 10” 3 
m/s 2 along x-axis and 6.50 x 10 -3 m/s 2 along both y and z-axis. Figure (5.77) shows 
the error in relative position, the maximum errors being 48.66 m along radial, 7.9 m 
along the track and 21.98 m along z-axis. Figure (5.78) shows the error in relative 
velocity. The maximum error in relative velocity in this case is 0.186 m/s along 
radial, 0.049 m/s along track and 0.097 m/s along the z-axis. 



Figure 5.76: Error in the Control Input: Case 16 
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Figure 5.77: Error in the Relative Position: Case 16 
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Case 17: t = 7r/2 n 


For specified time of t — n/2n = 1388.4s, figure (5.79) shows the error in control 
input trajectory. It can be seen that the maximum error in control input is 4.3 x 10” 3 
m/s 2 along x-axis and 5.25 x 10 -4 m/s 2 along both y and z-axis. Figure (5.80) shows 
the error in relative position, the maximum errors being 617.1 m along radial, 14.69 
m along the track and 59.39 m along z-axis. Figure (5.81) shows the error in relative 
velocity. The maximum error in relative velocity in this case is 0.48 m/s along radial, 
0.053 m/s along track and 0.171 m/s along the z-axis. 



Figure 5.79: Error in the Control Input: Case 17 
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Figure 5.80: Error in the Relative Position: Case 17 
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Case 18: t = n/n 


For specified time of t — n/n = 2776.8s, figure (5.82) shows the error in control input 
trajectory. It can be seen that the maximum error in control input is 6.49 x 10” 3 
m/s 2 along x-axis and 2.95 x 10 -3 m/s 2 along both y and z-axis. Figure (5.83) shows 
the error in relative position, the maximum errors being 740.6 m along radial, 310 m 
along the track and 845.2 m along z-axis. Figure (5.84) shows the error in relative 
velocity. The maximum error in relative velocity in this case is 1.21 m/s along radial, 
0.461 m/s along track and 0.97 m/s along the z-axis. 



Figure 5.82: Error in the Control Input: Case 18 
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5.8 Comparison of Optimal Control Input with 
and without the Effects of Earth’s Oblateness 

The following section deals with the cases when the effects due to oblateness of Earth 
are taken into account as discussed in Section 4.4. 

5.8.1 Varying the Initial Relative Velocity of the Target 

Case 19 

Rendezvous Problem: Consider a case, when a spacecraft is at a relative distance of 

173.2 Km from International Space Station (ISS) has to rendezvous with it within a 
period of (7r /n) seconds, where ISS is located at an approximated altitude of 400 
Km above Earth surface, having orbital frequency n = 1.1313658 x 10" 3 . 

We have the following data, 
x(t = 0) = aq(0) = 100 Km 
v x (t = 0) = £ 2 ( 0 ) = 0 Km/s 
y{t = 0) = £ 3 ( 0 ) = 100 Km 
v y (t = 0) = x 4 (0) = 0 Km/s 
z(t = 0) = x 5 (0) = 100 Km 
v z (t = 0) = x 6 (0) = 0 Km/s 
and 

tf = 7T fn s 

In this case, since, we are considering the effects of oblateness of Earth, it is required 
to take in account the inclination of the orbit of the target spacecraft, which in this 
case, is the ISS, having an inclination of 51.6° and the initial sum of the true anomaly 
and the argument of the perigee is 0°. 

The results of the simulation has been presented in the form of figures (5.85) to 
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(5.88). Figure (5.85) shows the propagation of error in the control input along the 
respective axes for the required control input, when the effects of Earth’s oblateness 
are introduced in the HCW model. Figure (5.86) plots the error in the relative 
distances along respective axes with time, whereas figure (5.87) plots the error in 
the relative velocities along the respective axes. Figure (5.88) shows the geometry of 
the maneuvering spacecraft orbit as predicted by the HCW model with and without 
oblateness effects, ft. can be observed that the basic shape of the orbit is same in 
both the cases and is in good agreement. 
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Figure 5.85: Error in the Control Input: Case 19 
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Figure 5.88: Geometry of the Maneuvering Spacecraft Orbit: Case 19 


Case 20 

Rendezvous Problem: Consider a case, when a spacecraft at a relative distance of 
173.2 Km from International Space Station (ISS) has to rendezvous with it within a 
period of (tt/ti) seconds, where ISS is located at an approximated altitude of 400 
Km above Earth surface, having orbital frequency n = 1.1313658 x 10" 3 . 

We have the following data, 
x(t = 0) = xi(0) = 100 Km 
v x (t = 0) = x 2 (0) = 0.05 Km/s 
y(t = 0) = x 3 (0) = 100 Km 
v y (t = 0) = x 4 (0) = 0.05 Km/s 
z(t = 0) = x 5 (0) = 100 Km 
v z (t = 0) = xe(O) = 0.05 Km/s 
and 

tf = 7T /n s 
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In this case, since, we are considering the effects of oblateness of Earth, it is required 
to take in account the inclination of the orbit of the target spacecraft, which in this 
case, is the ISS, having an inclination of 51.6° and the initial sum of the true anomaly 
and the argument of the perigee is 0°. 

The results of the simulation has been presented in the form of figures (5.89) to 
(5.92). Figure (5.89) shows the propagation of error in the control input along the 
respective axes for the required control input, when the effects of Earth’s oblateness 
are introduced in the HCW model. Figure (5.90) plots the error in the relative 
distances along respective axes with time, whereas figure (5.91) plots the error in 
the relative velocities along the respective axes. Figure (5.92) shows the geometry of 
the maneuvering spacecraft orbit as predicted by the HCW model with and without 
oblateness effects. It can be observed that the basic shape of the orbit is same in 
both the cases and is in good agreement. 



Figure 5.89: Error in the Control Input: Case 20 
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Figure 5.92: Geometry of the Maneuvering Spacecraft Orbit: Case 20 


Case 21 

Rendezvous Problem: Consider a case, when a spacecraft at a relative distance of 

173.2 Km from International Space Station (ISS) has to rendezvous with it within a 

period of (tt/ti) seconds, where ISS is located at an approximated altitude of 400 

Km above Earth surface, having orbital frequency n = 1.1313658 x 10" 3 . 

Considering only co-planar case, we have following data, 

x(t = 0) = xi(0) = 100 Km 

v x (t = 0) = x 2 (0) = 0.1 Knr/s 

y(t = 0) = x 3 (0) = 100 Km 

v y (t = 0) = x 4 (0) = 0.1 Km/s 

z(t = 0) = x 5 (0) = 100 Km 

v z (t = 0) = Xg(0) = 0.1 Km/s 

and 

tf — n/n s 
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In this case, since, we are considering the effects of oblateness of Earth, it is required 
to take in account the inclination of the orbit of the target spacecraft, which in this 
case, is the ISS, having an inclination of 51.6° and the initial sum of the true anomaly 
and the argument of the perigee is 0°. 

The results of the simulation has been presented in the form of figures (5.93) to 
(5.96). Figure (5.93) shows the propagation of error in the control input along the 
respective axes for the required control input, when the effects of Earth’s oblateness 
are introduced in the HCW model. Figure (5.94) plots the error in the relative 
distances along respective axes with time, whereas figure (5.95) plots the error in 
the relative velocities along the respective axes. Figure (5.96) shows the geometry of 
the maneuvering spacecraft orbit as predicted by the HCW model with and without 
oblateness effects. It can be observed that the basic shape of the orbit is same in 
both the cases and is in good agreement. 



Figure 5.93: Error in the Control Input: Case 21 
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Figure 5.96: Geometry of the Maneuvering Spacecraft Orbit: Case 21 


Case 22 

Rendezvous Problem: Consider a case, when a spacecraft at a radial distance of 173.2 
Km from International Space Station (ISS) has to rendezvous with it within a period 
of (vr/n) seconds, where ISS is located at an approximated altitude of 400 Km above 
Earth surface, having orbital frequency n = 1.1313658 x 10 -3 . 

We have the following data, 
x(t = 0) = xi(0) = 100 Km 
v x (t = 0) = £ 2 ( 0 ) = 0.2 Km/s 
y{t = 0) = £ 3 ( 0 ) = 100 Km 
v y (t = 0) = x 4 (0) = 0.2 Km/s 
z(t = 0) = x 5 (0) = 100 Km 
v z (t = 0) = X6(0) = 0.2 Km/s 
and 

tf = 7r/n s 
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In this case, since, we are considering the effects of oblateness of Earth, it is required 
to take in account the inclination of the orbit of the target spacecraft, which in this 
case, is the ISS, having an inclination of 51.6° and the initial sum of the true anomaly 
and the argument of the perigee is 0°. 

The results of the simulation has been presented in the form of figures (5.97) to 
(5.100). Figure (5.97) shows the propagation of error in the control input along the 
respective axes for the required control input, when the effects of Earth’s oblateness 
are introduced in the HCW model. Figure (5.98) plots the error in the relative 
distances along respective axes with time, whereas figure (5.99) plots the error in the 
relative velocities along the respective axes. Figure (5.100) shows the geometry of 
the maneuvering spacecraft orbit as predicted by the HCW model with and without 
oblateness effects. It can be observed that the basic shape of the orbit is same in 
both the cases and is in good agreement. 



Figure 5.97: Error in the Control Input: Case 22 
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Figure 5.100: Geometry of the Maneuvering Spacecraft Orbit: Case 22 


Table (5.6) summarizes the results of the cases (19-22). It can be observed that with 
the increase in the initial relative velocity of the maneuvering spacecraft, the errors 
in the control inputs, relative position and relative velocity as predicted by the HCW 
model, grows, when compared with the HCW model with the effects of oblateness of 
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the Earth. 
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5.8.2 Varying the Inclination of the Target Orbit 

Case 23 

Rendezvous Problem: Consider a case, when the maneuvering spacecraft is at a radial 
distance of 173.2 Km from the target which is in a circular orbit of radius 6778.14 
Km, and having an inclination of 0° and the initial sum of the true anomaly and the 
argument of the perigee is 0°. It is required to complete the low-thrust rendezvous 
in (7r jn) seconds, where n is the orbital frequency. 

We have following data, 
x(t = 0) = xi(0) = 100 Km 
v x (t = 0) = £ 2 ( 0 ) = 0.05 Km/s 
y(t = 0) = £ 3 ( 0 ) = 100 Km 
v y (t = 0) = £ 4 ( 0 ) = 0.05 Km/s 
z(t = 0) = £ 5 ( 0 ) = 100 Km 
v z (t = 0) = £6(0) = 0.05 Km/s 
and 

tf = 7 T /n s 

The results of the simulation has been presented in the form of figures (5.101) to 
(5.104). Figure (5.101) shows the propagation of error in the control input along the 
respective axes for the required control input, when the effects of Earth’s oblateness 
are introduced in the HCW model. Figure (5.102) plots the error in the relative 
distances along respective axes with time, whereas figure (5.103) plots the error in 
the relative velocities along the respective axes. Figure (5.104) shows the geometry of 
the maneuvering spacecraft orbit as predicted by the HCW model with and without 
oblateness effects. It can be observed that the basic shape of the orbit is same in 
both the cases and is in good agreement. 
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Figure 5.103: Error in the Relative Velocity: Case 23 



Figure 5.104: Geometry of the Maneuvering Spacecraft Orbit: Case 23 








150 


Case 24 

Rendezvous Problem: Consider a case, when the maneuvering spacecraft is at a radial 
distance of 173.2 Km from the target which is in a circular orbit of radius 6778.14 
Km, and having an inclination of 15° and the initial sum of the true anomaly and the 
argument of the perigee is 0°. It is required to complete the low-thrust rendezvous 
in (vr/n) seconds, where n is the orbital frequency. 

We have following data, 
x(t = 0) = aq(0) = 100 Km 
v x (t = 0) = £ 2 ( 0 ) = 0.05 Km/s 
y{t = 0) = £ 3 ( 0 ) = 100 Km 
v y (t = 0) = £ 4 ( 0 ) = 0.05 Km/s 
z(t = 0) = £ 5 ( 0 ) = 100 Km 
v z (t = 0) = x 6 (0) = 0.05 Km/s 
and 

tf = n/n s 

The results of the simulation has been presented in the form of figures (5.105) to 
(5.106). Figure (5.105) shows the propagation of error in the control input along the 
respective axes for the required control input, when the effects of Earth’s oblateness 
are introduced in the HCW model. Figure (5.106) plots the error in the relative 
distances along respective axes with time, whereas figure (5.107) plots the error in 
the relative velocities along the respective axes. Figure (5.108) shows the geometry of 
the maneuvering spacecraft orbit as predicted by the HCW model with and without 
oblateness effects. It can be observed that the basic shape of the orbit is same in 
both the cases and is in good agreement. 
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Figure 5.107: Error in the Relative Velocity: Case 24 



Figure 5.108: Geometry of the Maneuvering Spacecraft Orbit: Case 24 
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Case 25 

Rendezvous Problem: Consider a case, when the maneuvering spacecraft is at a radial 
distance of 173.2 Km from the target which is in a circular orbit of radius 6778.14 
Km, and having an inclination of 30° and the initial sum of the true anomaly and the 
argument of the perigee is 0°. It is required to complete the low-thrust rendezvous 
in (vr/n) seconds, where n is the orbital frequency. 

We have following data, 
x(t = 0) = aq(0) = 100 Km 
v x (t = 0) = £ 2 ( 0 ) = 0.05 Km/s 
y{t = 0) = £ 3 ( 0 ) = 100 Km 
v y (t = 0) = £ 4 ( 0 ) = 0.05 Km/s 
z(t = 0) = £ 5 ( 0 ) = 100 Km 
v z (t = 0) = x 6 (0) = 0.05 Km/s 
and 

tf = n/n s 

The results of the simulation has been presented in the form of figures (5.109) to 
(5.112). Figure (5.109) shows the propagation of error in the control input along the 
respective axes for the required control input, when the effects of Earth’s oblateness 
are introduced in the HCW model. Figure (5.110) plots the error in the relative 
distances along respective axes with time, whereas figure (5.111) plots the error in 
the relative velocities along the respective axes. Figure (5.112) shows the geometry of 
the maneuvering spacecraft orbit as predicted by the HCW model with and without 
oblateness effects. It can be observed that the basic shape of the orbit is same in 
both the cases and is in good agreement. 
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Figure 5.111: Error in the Relative Velocity: Case 25 



Figure 5.112: Geometry of the Maneuvering Spacecraft Orbit: Case 25 
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Case 26 

Rendezvous Problem: Consider a case, when the maneuvering spacecraft is at a radial 
distance of 173.2 Km from the target which is in a circular orbit of radius 6778.14 
Km, and having an inclination of 45° and the initial sum of the true anomaly and the 
argument of the perigee is 0°. It is required to complete the low-thrust rendezvous 
in (vr/n) seconds, where n is the orbital frequency. 

We have following data, 
x(t = 0) = aq(0) = 100 Km 
v x (t = 0) = £ 2 ( 0 ) = 0.05 Km/s 
y{t = 0) = £ 3 ( 0 ) = 100 Km 
v y (t = 0) = £ 4 ( 0 ) = 0.05 Km/s 
z(t = 0) = £ 5 ( 0 ) = 100 Km 
v z (t = 0) = x 6 (0) = 0.05 Km/s 
and 

tf = n/n s 

The results of the simulation has been presented in the form of figures (5.113) to 
(5.116). Figure (5.113) shows the propagation of error in the control input along the 
respective axes for the required control input, when the effects of Earth’s oblateness 
are introduced in the HCW model. Figure (5.114) plots the error in the relative 
distances along respective axes with time, whereas figure (5.115) plots the error in 
the relative velocities along the respective axes. Figure (5.116) shows the geometry of 
the maneuvering spacecraft orbit as predicted by the HCW model with and without 
oblateness effects. It can be observed that the basic shape of the orbit is same in 
both the cases and is in good agreement. 
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Figure 5.113: Error in the Control Input: Case 26 
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Figure 5.115: Error in the Relative Velocity: Case 26 



Figure 5.116: Geometry of the Maneuvering Spacecraft Orbit: Case 26 
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Case 27 

Rendezvous Problem: Consider a case, when the maneuvering spacecraft is at a radial 
distance of 173.2 Km from the target which is in a circular orbit of radius 6778.14 
Km, and having an inclination of 60° and the initial sum of the true anomaly and the 
argument of the perigee is 0°. It is required to complete the low-thrust rendezvous 
in (vr/n) seconds, where n is the orbital frequency. 

We have following data, 
x(t = 0) = aq(0) = 100 Km 
v x (t = 0) = £ 2 ( 0 ) = 0.05 Km/s 
y{t = 0) = £ 3 ( 0 ) = 100 Km 
v y (t = 0) = £ 4 ( 0 ) = 0.05 Km/s 
z(t = 0) = £ 5 ( 0 ) = 100 Km 
v z (t = 0) = x 6 (0) = 0.05 Km/s 
and 

tf = n/n s 

The results of the simulation has been presented in the form of figures (5.117) to 
(5.120). Figure (5.117) shows the propagation of error in the control input along the 
respective axes for the required control input, when the effects of Earth’s oblateness 
are introduced in the HCW model. Figure (5.118) plots the error in the relative 
distances along respective axes with time, whereas figure (5.119) plots the error in 
the relative velocities along the respective axes. Figure (5.120) shows the geometry of 
the maneuvering spacecraft orbit as predicted by the HCW model with and without 
oblateness effects. It can be observed that the basic shape of the orbit is same in 
both the cases and is in good agreement. 
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Figure 5.119: Error in the Relative Velocity: Case 27 



Figure 5.120: Geometry of the Maneuvering Spacecraft Orbit: Case 27 
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Case 28 

Rendezvous Problem: Consider a case, when the maneuvering spacecraft is at a radial 
distance of 173.2 Km from the target which is in a circular orbit of radius 6778.14 
Km, and having an inclination of 90° and the initial sum of the true anomaly and the 
argument of the perigee is 0°. It is required to complete the low-thrust rendezvous 
in (vr/n) seconds, where n is the orbital frequency. 

We have following data, 
x(t = 0) = aq(0) = 100 Km 
v x (t = 0) = £ 2 ( 0 ) = 0.05 Km/s 
y{t = 0) = £ 3 ( 0 ) = 100 Km 
v y (t = 0) = £ 4 ( 0 ) = 0.05 Km/s 
z(t = 0) = £ 5 ( 0 ) = 100 Km 
v z (t = 0) = x 6 (0) = 0.05 Km/s 
and 

t f = n/n 

The results of the simulation has been presented in the form of figures (5.121) to 
(5.124). Figure (5.121) shows the propagation of error in the control input along the 
respective axes for the required control input, when the effects of Earth’s oblateness 
are introduced in the HCW model. Figure (5.122) plots the error in the relative 
distances along respective axes with time, whereas figure (5.123) plots the error in 
the relative velocities along the respective axes. Figure (5.124) shows the geometry of 
the maneuvering spacecraft orbit as predicted by the HCW model with and without 
oblateness effects. It can be observed that the basic shape of the orbit is same in 
both the cases and is in good agreement. 
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Figure 5.121: Error in the Control Input: Case 28 
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Figure 5.123: Error in the Relative Velocity: Case 28 
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Figure 5.124: Geometry of the Maneuvering Spacecraft Orbit: Case 28 


It has been found that the error in the control inputs increases when the inclination 
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of the target orbit has been increased from 0° (equatorial orbit) to 15°, and so as the 
errors in the relative velocity and relative distance of the maneuvering spacecraft. 
However, when the inclination has been further increased to 30° and then to 45°, 
the respective errors decreases. As the orbital inclination increases to 60° and then 
to 90°, the errors in the control input, relative velocity and the relative distance of 
the maneuvering spacecraft are found to be decreasing. 

Table 5.7 summarizes the results of the cases (23-28). 
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Chapter 6 


Conclusions 


In this study, the unconstrianed optimal relative transfer in the orbit using low-thrust 
has been exhaustively studied by considering the model of Earth with approximated 
effects of oblateness and then, with the simplified model of spherical Earth. As 
discussed in the Chapter 1, the control input (specific thrust) for using low-thrust 
has to be as low as of the order of 10~ 4 Km/s -4 , and it can be seen from the results 
of the various numerical simulations presented in Chapter 5 that it indeed agrees 
with the requirement. 

Starting with the case of when the control input (Cases 1-2) is applied only along 
the in-track direction, it was observed that the requirements of the specific thrust 
were of the order of ICC 3 Km/s -2 , which is high as compared with the the low-thrust 
availability from the possible sources such as solar-sail. Hence, it is necessary to 
consider control inputs along both in-track and out-of the track, that is, coplanar 
inputs. 

In the case of the unconstrained coplanar optimal relative transfer (Case 3-5), it was 
observed that the required control input, that is specific thrust requirement comes 
down to the order of 10 -4 Km/s" 2 , which is well within the range of the low-thrust 
propulsion. 

To check the applicability of the approximated model [3] for the unconstrained ren- 
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dezvous trajectory optimization, the results of the optimization have been compared 
with those obtained by using the non-linear equations of relative motion in the orbit. 

The comparison has been done for the different cases of coplanar rendezvous in the 
orbit. In the optimal control problem, provided the other parameters constant and 
changing the initial velocity of the maneuvering spacecraft relative to the target- 
spacecraft from 0 Km/s to 0.2 Km/s, it has been found that the error in the predic¬ 
tion of the required control input trajectory increases as the initial relative velocity 
increases. In the same problem, provided the other parameters constant and changing 
the initial relative distance between the target and the maneuvering spacecraft from 
70.1 Km to 707.1 Km, it has been found that as the initial relative distance increases, 
the error in the required control input trajectory significantly increases. It has 
been observed provided the other parameters constant and changing the time of the 
rendezvous operation, the error in the required control input trajectory decreases as 
the time of the operation has been increases. The numerical results of the simulation 
have been tabulated in Table (5.1), (5.2), and (5.3). 

To check the applicability of the approximated model [3] for the unconstrained 
optimization of the station-keping problem, the results of the optimization have been 
compared with those obtained by using the non-linear equations of relative motion in 
the orbit. The comparison has been done for the different cases of optimal coplanar 
rendezvous in the orbit. From the simulated cases, it has been found that for different 
initial velocities of the maneuvering spacecraft relative to the target spacecraft, the 
required control input decreases with time. However, it is also observed that as the 
time of the operation has been increased the maneuvering spacecraft comes closer to 
the target than the specified nominal relative distance. Hence, to find the optimal 
time of operation for the optimal low-thrust unconstrained station-keeping mission, 
it is required to consider the restrictions on the time of operation imposed by the 
nominal relative distance requirements. Table (5.4) gives the optimal mission time 
for the station-keeping mission found in this way. 

To check the accuracy of the numerical solution which is obtained using MATLAB®function 
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bvp4c, with the closed-form solution of the optimization problem obtained in Chapter 
4, two cases have been compared. By superimposing the resulting plots of the optimal 
soltuion as obtained using the HCW on the plots obtained using the closed-form 
solution, it has been found that the results of the optimization performed numerically 
are in good agreement with the results obtained using the closed-form solution. 

To check the applicability of the approximate equations for the general relative 
motion in the orbit for the spherical Earth model, a number of cases have been 
studied for the same initial conditions and varying the time of the operation. It has 
been found that the error in the x-component of the control input increases as the 
time specified for the transfer increases, however, the error in the y- and z-component 
decreases significantly and then again increases. However, it can be observed that 
the errors in both the relative distance and the relative velocity increases with as 
the specified time of the optimal transfer increases. 

To account for the significant effects of the Earth’s oblateness, the perturbation 
effects due to J 2 have been added into the model. This presents an approximated 
oblate Earth model for the consideration. In the cases (19-22), the error that comes 
into the optimal control solution by introducing oblateness effect has been studied, 
by varying the initial relative velocity of the maneuvering spacecraft. For an specified 
inclination of the target orbit as 51.6°, which is the inclination of the ISS in the orbit, 
it has been found by varying the relative velocity of the maneuvering spacecraft from 
OKm/s 2 to 0.2Km/s 2 that the errors in the control input, in the relative velocity and 
in the relative distance, increases. 

Since, when considering the oblate Earth model, the inclination of the reference 
orbit is also a considerable parameter, in the cases (23-28), the numerical simulation 
has been carried out by varying the inclination of the reference orbit. It has been 
found by studying cases (23-28) that as the inclination of the reference orbit has 
been increased from 0° to 15°, the errors in the control input, relative velocity and 
the relative distance increases. However, the value of the control input decreases as 
the inclination of the reference orbit increase in steps from 15° to 90°, which denotes 
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the polar orbit. Tables (5.7) tabulates the respective results of the simulation found 
by including the approximate oblateness effect of Earth. 


6.1 Scope for Future Work 

This work demonstrates the capabilities of the low-thrust motion for the relative 
motion in the orbit. However, for more realistic model, the effects of oblateness has 
been considered but there are other significant effects which can be further added 
into the model, such as, the effects of atmospheric drag in the low-Earth orbits and 
gravitational effects of the Sun, moon and the other bodies in the space. Furthermore, 
in this study, starting with the Clohessy-Wiltshire model [3], and upt.o the oblate 
Earth relative motion model, the reference orbit was considered to be circular. The 
problem can be generalized by considering the optimization in case the reference 
orbit is elliptical, of which circular reference orbit can be a particular case. Since, the 
study has been concentrated on Earth, a more challenging problem can be using the 
similar approach for a heavenly body of arbitrary shape, with rather more complex 
zonal harmonics. 

This work presents the applicability of the approximated relative motion model 
for unconstrained low-thrust optimization and it can be further used to design a 
constrained optimal low-thrust control system. 
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Appendix A 


Gravitational Perturbations due 
to a Non-Spherical Planet 


The following derivation has been adapted from [14]. The gravitational force acting 
per unit mass on an object located at f relative to the center of the Earth is given 
by Newton’s law of gravitation as, 


/= (A.l) 

where, /i = 398600.4 Km 3 /s 2 , is the gravitational parameter of the planet Earth. 
Defining the potential function, 



r 


(A.2) 


such that, 


/ = V T 


(A.3) 
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where, 


_ „ <9 „ d „ <9 

v = x ^~ + 

ox ay oz 


The force per unit mass acting on a mass m 2 due to a series of point masses is given 
as, 


f- V Gm if 

J ~ / > r 3 


(A.4) 


which can likewise be obtained from the potential function as, 


X = E 


Grrij 


(A.5) 


for a continuous body, m, —* dm = pdV and yh —> 


v 


A.l Gravitational Potential for a Body of Arbi¬ 
trary Shape 


Consider an arbitrary shaped body of total mass m\. Let T be some reference frame 
with unit basis vectors x , y and z. Let dm(r > ) be a mass element in the body located 
at, 


r y = x'x + y'y + z' z (A.6) 

Using spherical coordinates, (V, A', 5'), 


x' = r' cos 5' cos 

(A.7) 

y = r cos 5' sin \' 

(A.8) 

z — r sin S 

(A.9) 



Ill 


To find the gravitational potential at a point r outside the body, with coordinates, 


r — xx + yy + zz (A. 10) 


with associated spherical coordinates, 


x — r cos S cos A 

(A.11) 

y = r cos S sin A 

(A.12) 

z = r sin S 

(A.13) 




Substituting equation (A.15) in equation (A.14), we get, 

dT( 

So that, the total gravitational potential becomes, 




(A.15) 


(A.16) 
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In which, 


1+7 — 2 ( 7 ) COS 7 


r ( r \ 

P 0 cosyj H-Pi cosyj + — P> cos7 + ... 

r \ r / 


OO / /\ n 

^( 7 ) ^(cos7) 

n=0 ' ^ 


(A. 18) 


where, P„(T) are the Legendre Polynomials stated as, 


(n + 1 )P n+1 (x) = (2 n + l)xP n (x) - nP n _i(x) 


and Pq(x) = 1, P\(x) = x, P 2 ^x) = |a; 2 — Therefore, 


OO n ' 

T(f) = — / P(0 ( — ) Pn(cOS^)dV 


-n -Z7 


n=0 


(A.19) 


Expanding upto n = 2, 


T(f) — — f p(r)dV — „ / p(r)r' cos^dV + — 


'7 


where, 


'7 


G r ij y P ^ (7) P n( COS ^) dV 


(A.20) 


p(r)dV = mi 


'7 


cos 7 = 






Thus, 


p{r)r' cos'ipdV = -[ p(r)r-rdV p(r')r dVr (A.21) 


'7 



V 


But, yy j v p(f)Pd,V defines the center of mass of the body. Hence, 


p(f)r' cos ifdV = 0 


(A.22) 


>v 


Hence, equation (A.20) becomes, 


^ Gm i G ^ 
T =-- + — 






_o JV 


n =2 


/\ n 


Pif) ( — ) P n (cOSlf)dV 


(A.23) 


Now, the term Pyf- is just the two-body potential for a point mass. Hence, the rest 
of the term denotes the potential resulting from the perturbative force per unit mass, 
that is, 


T p (r) = ? 


OO « 

e/m 

n=2 J y 


,\ n 


r y ) ( — ) P n (cosif)dV 


(A.24) 


Using the spherical coordinates, 


r ■ P = rP cos if 

= (r cos 8 cos Xx + r cos 8 sin A y + r sin 8z) ■ (P cos 8' cos X'x + P cos 8' sin X'y + P sin 8' 
= rr'(cos 8 cos 8' cos(A — A') + sin 8 sin 8') (A.25) 


so that, 


cos if = cos 8 cos 8' cos (A — A') + sin 8 sin 8' 


(A.26) 


Theorem A.l Addition Theorem for Spherical Harmonics, 


n ( — V 

P n (cos if) = P n ( sin 8)P n ( sin 5') + 2V ” — ! [A n ^ m Af n + B nm B' n ] (A.27) 

•' [n + m ! 

m=l v ' 
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where, 


A n , m = -P n ,m(sin 5) cos(mA) 

(A.28) 

A! n m = P n>m (sin S') cos(mA') 

(A.29) 

B n ,m = P„, m (sinA)sin(mA) 

(A.30) 

B n, m = Pn,m{sin5')sm(m\') 

(A.31) 

and the associated Legendre functions are P n , m (x) are, 


(jm 

Pn,m{x ) = (1 X 2 ) m/2 dxm P n (x) 

(A.32) 


Now, substituting equation (A.26) in equation (A.24) and making use of theorem 
(A.l), we get, 


T p (r) =- 


OO „ 

E »( 

,._o Jv 


/ / \ n 

^( r \ 

r — 


oo n / \ j 

P n (smS)P„(smS')dV + 2 " ~ - , 

( n + m ) ! 


X J ^ (7 J [ P n,m(AnS , )cos(mX')P njrn (smS)cos(mX) 
+ P n ,m (sin V) sin(mA , )P njm (sinA) sin(mA)] 


(A.33) 


Based on equation (A.33), the following factors are defined, 

J » = --ST— f P^)(r'T-P„(sm5')dV (A.34) 

rr e ' m l Jv 

C n ,m = --^r — 2 7“——77 [ p{P){r) n P n , m {s\Yi5')cos{m\')dV (A.35) 

n™mi (n + m)l J v 

Sn,m = --^r — 2 [ n , -T7 f p(P){r)P n ^ m (s'm d') sm(m\')dV (A.36) 

R£m i (-n + m)\J v 


where, R e is the normalizing radius of the rn \. In this case, it is the radius of the 
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Earth. Then, 


U(r) 


Grn i 


n =2 




P ni^ S ) 


Re 


^2 ^2 ( — ) P n ,m(sm8) [C n , m cos(mA) + S'™,™ sin(mA) 


(A.37) 


71=2 771=1 X 7 J 

where, the coefficients J n , C' n , m and S 1 ^,™ are determined experimentally from satellite 
observations. 

For rotationally symmetric bodies about z, we have, C n;rn = S n . rn = 0, and then, 


Cm / D \ n 

T p (r) = — ^2 Jn { — ) p n(sin8) (A.38) 

r n =2 V r / 

For planet Earth, the most dominant perturbing effect is the J 2 term effect, which is 
the result of Earth’s oblate shape. 

For planet Earth, 


J 2 = 1.083 x 10~ 3 


(A.39) 


So that, the perturbing potential including only the J 2 effects is given as, 


T 


p 





(A.40) 


A.2 Perturbation Force Per Unit Mass Due to J 2 


In ECI frame, the position vector of the spacecraft is given by, 


r = xx G + yi) G + zz G 


(A.41) 
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we have, 


sin 5 = - 
r 


(A.42) 


so that, equation (A.40) becomes, 


T 


p 



(A.43) 


Taking the gradient of equation (A.43), we have, 


d_ 

dx 



d_ 

dx 


(x 2 + y 2 + z 2 ) 3 / 2 


3 2x x 

2 ( x 2 + y 2 + z 2 ) b / 2 r 5 


(A.44) 


Similarly, 


d_ 

dy 




(A.45) 


d_ 

dz 




(A.46) 


and, 


d / z 2 \ d f z 2 \ 

dx\r 5 ) dx\(x 2 + y 2 + z 2 ) 5 / 2 ) 

2 ( 5\ 2x 

- Z V 2 ) {x 2 + y 2 + z 2 y/ 2 
xz 2 

= - 5 tv ( a - 47 ) 


Similarly, 


d_ 

dy 



-5 


yz 2 

/pT 


(A.48) 
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and, 


°(C)=2±- 5 * 


(A.49) 


Thus, the perturbation force is given as the gradient of the perturbing potential, 


F P = VT p 


— —/iJ2-Re 


15 xz 2 3x \ A 

YZ^ + 2^) Xa + 


+ 


15 z 3 z 3 z \ A 

YF + V + 2FJ ZG 


15 yz 2 3 y \ 

2 r 7 + 2r 5 ) VG 

(A.50) 


Recognizing that, r = xxq + yyc + and that, z — r ■ zq, equation (A.50) can be 
written as, 


f? _ 3 fiJ 2 ri 

~ 9 r.5 


5--^-1 If - 2(f • 2 g)z g 


In cylindrical coordinates, 


(A.51) 


r = rx o 


(A.52) 


and, 


zq = sin i sin 6 ■ xq + sin i cos 9 • i/o + cos i ■ zq 


(A.53) 


Hence, 


Zq ■ t = r sin i sin 6 


(A.54) 



X 


Substituting equations (A.53) and (A.54) in equation (A.51), we get, 


e? _3 yJ 2 R, 
p ~2 r 4 


(5 sin 2 i sin 2 9 — l) (rx 0 ) 


2r sin i sin 9 ( sin i sin 6 ■ x 0 + sin i cos 9 ■ y 0 + cos i ■ z 0 ) 


- _3^J 2 R 2 e 
p ~2 r 4 


(3 sin 2 i sin 2 9 — l) xq — (sin 2 i sin 29) y 0 — (sin 2 i sin 9) zq 


(A.55) 


(A.56) 



Appendix B 


Tensors of the Gravitational and 
Perturbation Forces 


B.l Tensor of the Perturbation Force 


The perturbation force can be represented in the (r, 6, i ) frame as, 


F p — (F r e r + Fgee + F^i) 

— — ^(T — 3sin 2 isin 2 d)e r + (sin 2 i sin 29)eg + (sin2f sind)e^ (B.l) 

The gradient of F p is a tensor, which can be represented as a dyadic product of the 
vector with the gradient operator as, 


^ d .Id . 1 d x 

v Fp yF r C r + F<)Cq + FiCi) (67*7— “h Cq —— + Ci 7 ~ 777 ) 

v or r t 06 r t sm 0 oi 


8F r 

1 dF r _ 

Fe 

1 dFr Ft 

dr 

rt. dO 

n 

rt sin 6 di rt 

dFg 

1 dFg , 

n dO " r 

Fr 

— — cot QF 

rt sm 6 oi rt 

dr 

n 

m 

dr 

i_dF\ 
rt dO 


—+ cot 9 Ft + 

rt sin 6 oi rt 


(B-2) 
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where, 


d F r 
dr 


6 J 2 /ii? 2 


(1 — 3 sin 2 i sin 2 8) 


(B.3) 


dF r 

~de 


9 J 2 h^e 

2 rf 


(sin 2 i sin 29) 


(B.4) 


d F r 


9 J 2 ^R 2 e 

2 rf 


(sin 2i sin 2 6) 


(B.5) 


Substituting the values of respective partial derivatives from equations (B.3), (B.4) 
and (B.5) in equation (B.2), we get, 


VF p = 


6 J 2 p-R, 


1 — 3 sin 2 i sin 2 8 
sin 2 i sin 26 

sin 2 i sin 8 


sin 2 i sin 26 

— \ + sin 2 if ^ sin 2 8 - \ 


— \ sin 2% cos 8 


sin 2 i sin 6 


— | sin 2 % cos 8 


■f + sin 2 i I | sin 2 6 + \ 


(B.6) 
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B.2 Tensor of the Gravitational Force 



g(r) = 4 e r (B.8) 

n 

So that, equation (B.7) becomes, 

0 0 

0 (B.9) 

0 




Appendix C 


Detailed Integration Steps 


The following integration steps were used in integrating the respective elements of 
equation (4.67). 


sin(nr) sin(n(r — t))dr = j ^ — - cos(2nr — nt) + - cos(nf)^ dr 


=-sin(2nr — nt) 4—r cos (nt) 

4 n 2 


(C.l) 


'1 1 

cos (nr) sin(n(r — t))dr = / ( - sin(2nr — nt) + - sin (nt) ) dr 


>)■ 


1 1 

=-cos(2 nr — nt) -r sin(nf) (C.2) 

4n 2 


1 1 

sin(nr) cos(n(r — t))dr = \ ( - sin(2nr — nt) + - sin (nt) ) dr 


>> 


=-cos(2nr — nf) 4—r sin(nf) 

4n 2 


(C.3) 
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cos(nr) cos(n(r — t))dr — j cos(2nr — nt) + - cos(nt)^ dr 


1 , 1 , , 

= — sin(2nr — nt) H—rcos(nf) 
An 2 


(C.4) 


If 1 

rcos(n(r — t))dr = -/ sin(n(r — t))dr A —rsin(n(r — t )) 

Tl J 71 

1 1 

= — cos(n(r — t)) H—rsin(n(r — t )) 


rr 


n 


(C.5) 


From Eq. (4.41), we get, 


66 + -66 dr 
n 


—^ ( sin(n(f — r)) (2C471 sin(nr) — Ci sin(nr) + C 2 n cos(nr) + 2C3 cos (nr) — 2C3) 


-2 ( (2 cos (n(£ — r)) — l) (2C 2 n sin(nr) + 4C3 sin(nr) — 4C4U, cos (nr) + 2Ci cos (nr) 


rr 


3 Csnr + 3C471 — 2Ci)^ 


dr 


(C.6) 


where, 


j sin(n(t — r)) {2C^n sin(nr) — C\ sin(nr) + C 2 n cos (nr) + 2C 3 cos(nr) — 2 C 3 )dr 
= (2 C 4 n — Ci) [ sin(nr) sin(n(r — f))dr + (C 2 n + 2C 3 ) [ cos (nr) sin(n(r — t))dr 


2C3 / sin(n(r — f))dr 


(C.7) 


Substituting values of the respective integrals from equations (C.l) and (C.2) in 
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(C.T) and evaluating the rest of the integrals, we get, 


sin (n(t — r)) (2 C 4 n sin(nr) — C\ sin (nr) + C 2 n cos (nr) + 2 C 3 cos (nr) — 2 C 3 )dr 


1 T 

= (2nC 4 — Ci) ( — — sin(2nr — nt) H— cos (nt) 
4 n 2 

+ (C 2 n + 2 C 3 ) ( — — cos(2nr — nt) — -r sin(nt)^ + 2 — cos(n(r — t)) 


4 n 


n 


(C. 8 ) 


and, 


( 2 C 2 n + 4C 3 ) sin(nr) + (2Ci — 4C 4 n) cos (nr) — 3C 3 nr + 3C 4 n — 2C 4 


(cos(n(r — t)) — l)dr 

= 2 (C 2 n + 2 C 3 ) J sin(nr) cos(n(r — t))dr — 2 ( 2 C 4 n — Ci) J cos (nr) cos (n(r — t))dr 

— 3C 3 n J r cos(n(r — t))dr + (3C 4 n — 2 Ci) J cos (n(r — t))dr 

— 2 (C 2 n + 2C 3 ) f sin (nr) dr + 2(2C 4 n — C — 1) [ cos (nr) dr 


+ 3C 3 n J rdr + ( 2 C 4 — 3C 4 n) J dr 

— -^-( 2 C 4 n — Ci) sin(n( 2 r — t)) — -^-(C 2 n + 2C 3 ) cos(n(2r — t)) 

1 L Zj /1 

— 3C 3 rsin(n(r — £)) H—(3C 4 n — 2Ci) sin(n(r — £)) — — C 3 cos(n(r — t)) 

n n 

H—(2C 4 n — Ci) sin(nr) H—(C 2 n + 2C 3 ) cos (nr) H—C 3 nr 2 
n n 2 

+ (C 2 n + 2C 3 )rsin(nt) — (2C 4 n — Ci)rcos(nr) + (2Ci — 3C 4 n)r 


(C.9) 
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Substituting values from Eqs. (C. 8 ) and (C.9) into Eq. (C. 6 ), we obtain: 

J (66 + ^66^r 

= -^t( 2 C 4 n - Ci) sin(n( 2 r — t)) + -^(C 2 n + 2C 3 ) cos(n(2r - t)) 

4 n 6 4 n 6 

+ -^rC 3 r sin(n(r — t)) --(3C 4 n — 2 C 4 ) sin (n(r — t)) 

n 2 n 6 

6 2 

H—o C 3 cos(n(r - f)) H—-C 3 cos(n(r - t)) 
n 6 n 6 

--(2 C 4 n — Ci) sin (nr) --(C 2 n + 2C 3 ) cos (nr) — — C 3 r 2 

n d rr n 

5 5 

- ^ 2 (C 2 n + 2C 3 )rsin(ni) + ^^( 2 C 4 n - Ci)rcos(nt) 

—^-t( 2 Ci — 3C 4 n) + const. (C. 10 ) 

n 2 


From Eq.(4.42), we get, 


66 - 266 


^(4C 2 n + 8 C 3 ) sin(nr) + (4 Cl — 8 C 4 n) cos (nr) 


— 6 C 3 nr — 6 C 4 n — 4Ci ) sin(n(r — £)) + ( (2C 4 n — Ci) sin(nr) 


+ (C 2 n + 2C 3 ) cos (nr) — 2C — 3 ) cos(n(r — t)) 


(C.ll) 


and, 


(66 - 266) dr 


n 


(C 3 (4 sin (nr) — 3 nr) + 2C 2 nsin(nr) + 2Ci(cos(nr) — 1) 


+ C 4 n(3 — Acosinr))) sin(n(r — t)) ) dr — 

/ n 


(Ci — 2C 4 n) sin(nr) 

— C 2 ncos(nr) — 2C — 3(cos(nr) — 1) ) cos(n(r — t))dr (C. 12 ) 
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Now solving, 


(C 3 (4sin(nr) — 3 nr) + 2C' 2 nsin(nr) + 26 i(cos(t 7 t) — 1 ) 


+ 6477(3 — 4 cos(nr))) sin(n(r — t))J dr 
- AC^n 2 sin(nr) + 2C'insin(nr) 


1 

n 


+ 6 3 (3t7 — Ancos(nr )) — 3C 2 n 2 cos (nr)) cos (n(r — £))dr 
— — (C 3 (4sin(nr) — 3nr) + 2C 2 nsin(nr) + 26 i 77 (cos(t 7 t) — 1) 
+ 6 * 471(3 — 4cos(t7t))) cos(n(r — t)) 


(C.13) 


where, 


1 

n 


AC^n 2 sin(nr) + 2C\n sin(nr) 


+ 6 3 (3t 7 — 4ncos(nr)) — 3 C 2 n 2 cos(nr)) cos(n(r — t))dr 
-(2C\ — 46477) J sin( 77 r) cos(n(r — t))dr 

— (26 2 77 + 46 3 ) J cos (nr) cos(t7(t — t))hr + 36 3 J cos (n(r — t))dr 

- - ((2C 2 n — 46 3 ) sin(277T — nt)) H -((26i — 46477 ) cos(ti( 2 t — £))) 

477 477 

— —6 3 sin(77(r — t)) — — (6 3 (4sin(77r) — 3nr) 

n n v 

+ 26 2 77sin(77r) + 26 i(cos(t 7 t) — 1 ) + 6471(3 — 4 cos(t7t))) cos(ti(t — t)) 

— -r(26i — 46 4 77) sin(77t) + -r(26 2 77 — 46 3 ) cos(nt ) 


(C.14) 
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(Ci — 2C 4 n) J sin(nr) cos(n(r — t))dr — (C 2 + 2C 3 ) J cos(nr) cos(n(r — t))dr 
4-2C 3 J cos (n(T — t))d,T 

- — -'-((2C' 3 + C 2 n) sin(n(2r — £))) — ^-((C 4 — 2C 4 n) cos(n(2r — £))) 

2C 1 1 

4-- sin(n(r — t )) + -r(Ci — 2C — 4n) sin(nt) — -r(C 2 n + 2C 3 ) cos(nt) 

77 / ^ ^ 

(C.15) 


Substituting Eqs. (C.13)-(C.15) into Eq. (C.12), we get, 


(6<£e + 2fi &)dr 

— — — ((2C 2 n — 4C 3 ) sin(2nr — nt)) 

on /I o') ' ' 


n 4n v ' v " 

+ — ((2Ci — 4C 4 n) cos(n(2r — t))) — —C 3 sin(n(r — t)) 

-(C 3 (4sin(nr) — 3 nr) + 2C 2 nsin(nr) + 2Ci(cos(nr) — 1) 

X 

+ C 4 n(3 — 4 cos (nr))) cos (n(r — t )) — -r(2Ci — 4C 4 n) sin(nt) 

4—r(2C 2 n — 4C 3 ) cos (nt)-(C 3 ( 4 sin(nr) — 3nr) + 

2 n v 

2C 2 nsin(nr) + 2Cin(cos(nr) — 1) + C 4 n(3 — 4cos(nr))) cos(n(r — t)) + 

^ _ ^“(( 2C 3 + C2 ‘ n ) sin ( n ( 2r ~ *))) 

1 2(7 

-((Ci — 2C 4 n) cos(n(2r — t))) 4-- sin(n(r — t)) 

4 n n 

4- —t(Ci — 2C — 4n) sin (nt) — -r(C 2 n 4- 2C 3 ) cos (nt) (C 


(C.16) 



XX 


From Eq.(4.43), we get, 


- fi (3T -3 1 + ^ 3 ) ^ dr 

J ^3r H — sin (n(t — r)) — 3t^ ^C 3 ( — — sin (nr) + 3r) 

2 2 \ 

2 C 2 sin (nr) + Ci (-cos(nr)) + C 4 (4 cos (nr) — 3) I dr 

n n / 


2 2 r ^^ 

-cos nr — sin nr 

n n / \ n 


.22 . \ 

— 2C 4 sin(nr) + C 3 (-cos (nr)) — C 2 cos (nr) dr 

n n / 

= — — J (C 3 (4sin(nr) — 3nr) + 2C 2 nsin(nr) + 2Ci(cos(nr) — 1) 

+ C 4 n(3 — 4 cos(nr))) (4sin(n(r — t)) — 3nr + 3 nt)dr 

+ — J ((Ci — 2C 4 n) sin(nr) — C 2 ncos(nr) — 2C 3 (cos(nr) — 1)) ( cos(n(r — t)) — l)dr 

(C.17) 


in which, we have, 


J (C 3 (4sin(nr) — 3nr) + 2C 2 nsin(nr) + 2Ci(cos(nr) — 1) 

+ C 4 n(3 — 4 cos(nr))) (4 sin(n(r — t)) — 3nr + 3nt)dr 
=4C 3 /(4sin(nr)-3nr)sin(n(r-t))dr + 8C 2 n J sin(nr) sin(n(r — t))dr 

— 4C 4 n J (4 cos (nr) — 3) sin(n(r — t)) + 8C1 J (cos (nr) — 1) sin(n(r — t))dr 

— 3C 3 n J r(4sin(nr) — 3 nr)dr — QC-^n 2 J rsin(nr)dr 

+ (6C 2 n 2 t + 12C 3 ni) [ sin(nr)dr + 3C 4 n 2 [ r(4 cos(nr) — 3)dr 


+ (6C 4 nt — 12C 4 n 2 t) J cos(nr)dr — 6C 4 n J r(cos(nr) — l)dr 
— 9C 3 n 2 t [ rdr + (9C 4 n 2 t — 6C 4 nt) [ dr 


(C.18) 
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where, 


and, 


and, 


J (4sin(nr) — 3nr) sin(n(r — t))dr 

r i X 

= — / — (3n — An cos(nr)) cos(n(r — t))dr -(4sin(nr) 

J n n 

— 3 nr) cos (n(r — t)) 

1 3 

= — sin(n(2r — t)) -sin(n(r — t)) + 2 r cos (nt) 

n n 

X 

-(4sin(nr) — 3 nr) cos(n(r — t)) 

n 


J (4cos(nr) — 3) sin(n(r — t))dr 

=4 J cos (nr) sin(n(r — t))dr — 3 J sin(n(r — t))dr 
1 3 

=-cos(n(2r — t)) -\— cos (n{r — t)) — 2 r sin (nt) 

n n 


J (cos (nr) — 1) sin(n(r — t))dr 

~ J o° s ( nr ) s * n ( n ( T — t))d-T — J sin(ra(r — t))dr 

1 11 

=-cos(n(2r — t)) 4— cos(n(r — t)) -r sin (nt) 

4 n n 2 


and, 


J r(4sin(nt ) — 3 nr)dr 
=4 J r s'm(nr)dr — 3n J r 2 dr 

X 1 

=4(—sin(nr)-rcos(nr)) —nr 3 

\ r) ' 


n 


(C.19) 


(C.20) 


(C.21) 


(C.22) 



Substituting Eqs. (C.19)-(C.22) into Eq. (C.18), we obtain: 


J (C 3 ( 4 sin(nr) — 3 nr) + 2 C 2 nsin(nr) + 2 Ci(cos(nr) — 1 ) 

+ C 4 n(3 — 4 cos(nr))) (4 sin(n(r — t )) — 3nr + 3 nt)dr 

4 2 

=—C 3 sin(n( 2 r — t )) — 2 C 2 sin(n(2r — t ))- C\ cos(n( 2 r — t )) 

n n 

12 4 

+ 4C 4 cos(n(2r — t ))- C 3 sin (n(r — £)) H—C 3 (4sin(nr) — 3nr) cos(n(r — t )) 

n n 

g 

H— C\ cos(n(r — £)) — 12C 4 cos(n(r — £)) + 12C 4 nr sin(nr) 

12 

— 6 Cirsin(nr) + 6 Cit — 12C 4 nt)- C 3 sm(nr) — 6C 2 sin (nr) + 6 C 2 nrcos(nr) 

n 

6 

+ I 2 C 3 T cos (nr) — (6 C 2 nt + I 2 C 3 C cos (nr)-Ci cos (nr) + 12C 4 cos (nr) 

n 

+ 3 C 3 n 2 T 2 — ^C 3 n 2 tr 2 — ^C 4 n 2 r 2 + 3Ci?rr 2 + 8 C 4 nrsin(nt) 

— 4CiTsin(nt) + 4rC 2 n cos(nt) + 8C 3 t cos(nt) + ( QC^nH — 6Cint)r (C.23) 


And, 


j ((Ci — 2 C 4 / 7 ,) sin(nr) — C 2 n cos(nr) — 2C 3 (cos(nr) — 1)) ( cos(n(r — £)) — l)dr 
-{Ci — 2 C 4 n) f sin(nr) cos(n(r — t))dr — ( C 2 n — 2C 3 ) [ cos (nr) cos(n(r — t))dr 


+ 2 C 3 y cos(n(r — t))dr + ( 2 C 4 n — Ci) J sin(nr)dr + (C 2 n + 2 C 3 ) J cos(nr)dr 

= —(C 2 n + 2 C 3 ) sin(n(2r — £)) — — (Ci — 2 C 477 ,) cos(n( 2 r — t))+ 

—C 3 sin(n(r — £)) H—(C 2 n + 2C 3 ) sin(nr)-( 2 C 4 n — Ci) cos(nr) 

n n n 

+ -(Ci — 2C 4 n)r sin(nt) — -(C 2 n 2C 3 )t cos(nt) — 2 C 3 t 


(C.24) 
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Substituting Eqs.(C.23) and (C.24) in Eq.(C.17). we obtain: 


66 — 6 (3t — 3 1 -\ 

v n 


dr 


= - -^(C 2 n + 2 C 3 ) sin(n( 2 r - t )) + ~^C 2 sin(n( 2 r - t )) 

2 rr rr 

- -^C 3 sin(n(2r — t)) — -^(Ci - 2C 4 n) cos(n(2r - £)) 

77/ Z72 

4 2 12 

- 2 C 4 cos(n(2r — £)) 4—-Ci cos(n( 2 r — £)) 4—-C 3 sin(n(r — t )) 


rr 


rr 


rr 


4 12 

4—-C 3 (4 sin(nr) — 3nr) cos(n(r — £)) 4—-C 4 cos(n(r — t )) 


n° 


rr 


4 12 

—Ci cos(n(r — t)) 4 —-C 3 sin(n(r — t ))-C 4 r sin(nr) 


rr 


rr 


n 


+ -^-CiT sin(nr) — (-^-CC-C 4 t) sin(nr) 

rr v rr n ' 

2 6 12 
4 —^(2C 3 + C 2 n ) sin(nr) 4 —-C 2 sin(nr) 4 —-d 3 sin(nr) 


rr 


rr 


rr 


— — C 2 t cos (nr) --rC 3 cos(nr) + (— C 2 t 4 --C 3 t) cos (nr) 

n rr v n rr 7 

2 12 
--(2C 4 n — Ci) cos(nr)--C 4 cos (nr) 


rr 


rr 


+ —77C1 cos(nr) - 3 C 3 r 3 4 - |c 3 tr 2 
rr 2 

— — Cir 2 4- -C 4 r 2 4- -)t(Ci — 2 C 4 n)r sin(rrt) 
n 2 rr 

8 4 1 

-C 4 rsin(nt) 4 —-CiTsin(ni)--(C 2 n 4 - 2 C 3 )rcos(nt) 

n rr rr 

- C 2 t cos(nt) --C 3 T cos(nt) — (9C 4 t-Ci£)r 


n rr 

4 

—C 3 r + const, 
rr 


n 


(C.25) 


From Eq.(4.44), we get, 


J (-266-6(46 -3))dr 

- ^ ((Ci — 2C 4 n) sin(rrr) — C 2 n cos(nr) — 2C 3 (cos(nr) — 1)) sin(rr(r — t))dr 
' 3 ( 4 sin(nr) — 3nr) 4 - 2C 2 nsin(rrr) 4 - 2Ci(cos(nr) — 1) 

+ C 4 n(3 — 4cos(nr)))(4cos(n(r — t )) — 3)dr (C.26) 


1 

n 
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where, 


J ((Ci — 2C 4 n) sin(nr) — C 2 ncos(nr) — 2C 3 (cos(w) — 1)) sin(n(r — t))dr 
-(Ci — 2 C 4 n) f sin(nr) sin(n(r — t))dr + (—C 2 n — 2C 3 ) f cos (nr) sin(n(r — t))dr 


+ 2C 3 J sin (n(r — t))dT 

— —(Ci — 2C 4 n) sin(n(2r — £)) + ^-(2C 2 n + 2C 3 ) cos(n(2r — £)) 

2 1 1 

-C 3 cos(n(r — t)) + -(C 2 n + 2C 3 )rsin(ni) + -(C 4 — 2C 4 n)r cos(ni) 

72 ^ ^ 


(C.27) 


and, 


J (C 3 (4sin(nr) — 3nr) + 2C 2 nsin(nr) + 2Ci(cos(nr) — 1) 

+ C 4 n(3 — 4cos(nr)))(4cos(n(r — £)) — 3)dr 
=4C 3 J (4 sin(nr) — 3nr) cos(n(r — t))dr + 8C 2 n J sin(nr) cos (n(r — t))dr 

— 4 C 4 n J (4cos(nr) — 3) cos(n(r — t))dr + 8Ci J (cos (nr) — 1) cos (n(r — t))dr 

— (6 C 2 n + 12C 3 ) f sin (nr) dr + (12C 4 n — 6Ci) [ cos (nr) dr + 9C 3 n [ rdr 


+ (9Ci — 9C 4 n) / dr 


(C.28) 


and, 


(4sin(nr) — 3 nr) cos(n(r — £))dr 


=4 J sm(nr) cos(n(r — t))dr — 3n J rcos(n(r — t))dr 

1 1 

= — cos(n(2r — t)) H—(4sin(nr) — 3 nr) sin(n(r — £)) 
n n 

3 

— — cos(n(r — t)) + 2r sin(nt) 


(C.29) 
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Substituting Eqs. (C.27)-(C.29) into Eq. (C.28), we obtain: 


J ( 63(4 sin (nr) — 3 nr) + 26* 2 n sin(nr) + 2 C 1 (cos(nr) — 1 ) 

+ 6* 4 n(3 — 4cos(nr)))(4cos(n(r — t)) — 3 )dr 

2 4 

=—Ci sin(n( 2 r — t)) — 4 6* 4 sin(n( 2 r — t)) H— 6*3 cos(n( 2 r — t)) 
n n 

4 

— 2C 2 cos(n( 2 r — t)) H—C , 3 ( 4 sin(nr) — 3 nr) sin(n(r — t)) 

8 12 
- Ci sin(n(r — t)) + 126* 4 sin(n(r — t)) H - 6*3 cos (n(r — t)) 


n n 

H—(126' 4 n — 6 C 1 ) sin(rir) H—( 66 * 2 n + 126* 3 ) cos (nr) H— C 3 nr 2 
n n 2 

+ 4C 2 nr sin (nt) + 8 x 6*3 sin (nt) — 8 6 * 4 nr cos (nt) + 46*ir cos(nt) + r( 66 *i — 


96' 4 rt) 

(C.30) 


Substituting Eqs. (C.27)-(C.30) into Eq. (C.26), we obtain: 

J (-266-6(4^-3))dr 

1 4 

= — —— (C*i — 26* 4 n) sin(n( 2 r — t)) - 6 * 4 sin(n( 2 r — t)) 

ZjTX Tl 

+ -rCi sin(n( 2 r - £)) + 2 — ( 6 * 2 n + 26* 3 ) cos(n( 2 r - £)) 
rr 2 rr 

2 4 4 

- 6* 2 cos(n( 2 r — t)) H—- 6*3 cos(n( 2 r — £)) H—- 6 * 3 ( 4 sin(nT) — 3nr) sin(n(r — £)) 

n n z n z 

H -6* 4 sin (n(r — £)) —— 6* 4 sin(n(r — £))--6*3 cos(n(r — £)) 

n rr n z 

4 1 1 

--6*3 cos(n(r — t)) H- -(126* 4 n — 66 *i) sin(nr) H— -( 66 * 2 n + 126 * 3 ) cos (nr) 

n z rr n z 

9 1 8 

H— 6 * 3 r 2 H—( 6 * 2 n + 26* 3 )r sin(ni) H— 6 * 3 r sin(ni) + 46* 2 r sin(ni) 

2 n n 

1 4 

H—( 6 *i — 26* 4 n)r cos(nt) H— 6 *it cos(nt) — 86 * 4 rcos(ni) 
n 

H—( 66 *i — 96* 4 n)r + const. 
n 


(C.31) 
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The integral of Eq. (4.45) is obtained as: 


- / 6i£ 4 dr 


n 


If 1 

=-/ sin (nit — r)) (C 6 cos (nt) - C 5 sin (nr)) dr 

n J v n 

=- Cq f sin(n(t — r)) cos(nr)dr H— -=C 5 f sin(nr) sin(n(t — r))dr 


n 


rr 


4n 3 


n( — 2 r(C 6 nsin(nt) + C 5 cos (nt)) — Cq cos(n( 2 r — t )) 


+ —-C 5 sin(n( 2 r — t)) + const. 


The integral of Eq. (4.46) is obtained as: 

- J )dr 

1 r 1 

= — / cos (n(t — t)) (Cq cos (nr)-C 5 sin(nr))dr 

n J n 

=— C${-t sin(nt)-cos(n( 2 r — £))) 

n 2 4n 

1 1 

— Cq( —— sin(n( 2 r — t)) H—rcos(nt)) + const. 
v 4 n 2 


(C.32) 


(C.33) 
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